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ABSTRACT 


The  radiation  transport  subroutines  of  the  SPUTTER  code  for  spherical  geometry 
have  been  revised.  The  DIFFU  subroutine  has  been  eliminated  and  RADTN,  TRANS, 
and  KAPPA  have  been  recoded.  The  results  of  this  work  on  the  codes  are  (1) 
improved  logical  organization,  (2)  more  efficient  and  rapid  calculation,  (3) 
improved  accuracy,  (4)  more  complete  documentation,  and  (5)  comparisons  with 
test  problems.  The  much  simpler  and  more  accurate  diffusion  approximation  is 
exploited  when  a  new  diffusion  criterion  is  satisfied  in  shells  within  the 
sphere.  A  more  accurate  angular  integration  of  the  intensity  makes  use  of  the 
y-line  integration  results  more  efficiently  to  give  improved  fluxes.  Reorgani¬ 
zation  of  the  calculation,  saving  of  quantities  to  be  used  again,  and  use  of 
a  fast  exponential  routine  have  resulted  in  speeding  up  the  routines  by  approxi¬ 
mately  a  factor  of  2.  The  diffusion  and  angular  integration  improvements 
apparently  have  resulted  in  an  additional  factor  of  2  speedup  for  the  same 
accuracy. 
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The  SPUTTER  code  subroutines  for  radiation  transport  in  spheres  described 
herein  are  as  they  existed  on  July  1,  1965.  General  Atonic  has  exercised  due 
care  in  preparation,  but  does  not  warrant  the  merchantability,  accuracy,  and 
completeness  of  these  subroutines  or  of  their  description  contained  herein. 

Hie  complexity  of  this  kind  of  program  precludes  any  guarantee  to  that  effect. 
Therefore,  any  user  must  make  his  own  determination  of  the  suitability  of  these 
subroutines  for  any  specific  use  and  of  the  validity  of  the  information  produced 
by  their  use. 
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SECTION  I 


INTRODUCTION 


The  spherical  version  of  the  radiation  transport  routines  for  SPUTTER^*) 
has  been  completely  revised  to  increase  both  the  speed  and  the  accuracy  of 
the  calculation.  The  new  subroutines  contain  a  revised  diffusion  criterion, 
the  use  of  Planck  means  in  the  transparent,  or  thin,  limit  (Rosseland  means 
have  always  been  used  in  the  thick  limit),  and  an  improved  method  of 
integrating  over  angles. 

The  logic  in  the  SRADTN  subroutine  has  been  reorganized:  first,  to 
select  an  initial  set  of  y-lines  from  which  the  set  for  each  frequency  group 
will  be  chosen  later;  second,  to  merge  frequency  groups  for  which  the  source 
func  .ion  is  small,  even  for  the  highest  temperature  in  the  mesh;  third,  to 
limit  the  source  calculation  to  regions  of  appreciable  temperature  and  thereby 
reduce  the  transport  calculation;  and,  finally,  to  establish  diffusion  regions 
and  perform  the  diffusion  calculation  in  SRADTN  (formerly  this  was  done  in 
a  subroutine  called  DIFFU).  In  essence,  therefore,  the  new  SRADTN  sub¬ 
routine  establishes  the  limiting  set  of  active  zones  for  sources,  y-lines,  and 
frequency  groups.  The  new  STRANS  subroutine  contains  the  intensity  inte¬ 
gration  along  y-lines  for  the  frequency  groups  selected  in  SRADTN, 

What  has  been  achieved  in  receding  these  subroutines  has  been  the 
removal  from  inside  the  main  y-line  loop  the  redundant  calculations  made 
previously,  and  this  was  principally  done  by  storing  the  complete  set  of 
x-y  intersections.  The  improvement  in  the  angular  integration  is  achieved 
by  storing  intensities  along  the  previous  y-iine  and  interpolating  between 
y-lines  for  an  improved  quadrature  summation  of  fluxes.  Both  the  quadratic 
and  linear  forms  of  interpolation  were  tested,  but  since  the  linear  form 
appears  to  be  more  appropriate  in  most  cases,  it  is  used  in  STRANS. 

Appropriate  expansions  in  optical  depth  for  transparent  regions  as 
well  as  diffusion  boundary  conditions  are  used  in  the  new  subroutines.  At 
present  the  Rosseland  mean  is  used  in  the  thick  limit  and  the  Planck  mean 
in  the  thin  limit.  This  approach  is  only  temporary,  as  eventually  the 
more  realistic  transmission  means,  which  will  automatically  limit  correctly, 
will  be  used. 

The  numerical  solution  of  the  radiation  transport  equation  along 
selected  sampling  rays  through  a  sphere  is  discussed  in  Section  II.  The 
improved  diffusion  approximation  is  described  in  Section  III,  and  the 
frequency  integrations  needed  for  the  SPUTTER  calculations  are  described 
in  Section  IV. 
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Although  most  of  this  report  is  concerned  with  the  improved  SRADTN 
and  STRANS  routines,  some  changes  have  also  been  made  in  certain  auxilliary 
subroutines,  e.  g. ,  a  faster  exponential  subroutine  FREXP,  the  Planck 
function  subroutine  PLNKUT,  the  absorption  coefficient  interpolation  sub¬ 
routine  KAPPA,  etc.  These  improvements  are  discussed  in  Section  V,  and 
the  listings  of  the  subroutines  SRADTN,  STRANS,  FREXP,  KAPPA,  and 
PLNKUT  are  given  in  the  Appendixes.  The  preliminary  results  of  the  improve¬ 
ments  in  accuracy  and  calculation  times  from  comparison  with  the  old  sub¬ 
routines  are  given  in  Section  VI. 

The  obvious  differences  between  the  new  and  the  old  subroutines  are 
in  the  new  angular  integration  scheme  and  criteria  for  diffusion.  The  y-line 
integration  of  intensities  is  essentially  the  same.  Actually,  what  has  been 
achieved  is  a  careful  re-evaluation  of  the  complete  code,  which  has  resulted 
in  many  time  savings  and  in  some  increased  accuracy.  Additional  accuracy 
can  be  achieved  by  going  to  more  complete  and  therefore  more  complex 
integration  schemes  requiring  additional  core  storage.  These  schemes 
will  probably  need  the  increased  fast- storage  capacity  of  the  new  generation 
of  computers  (e.  g. ,  CDC  6600  and  IBM  360).  Careful  comparisons  have 
been  made  between  single  cycles  of  multifrequency  transport  carried  out 
with  the  old  code  and  with  the  new  code  (as  discussed  in  Section  VI).  These 
calculations  have  pointed  up  certain  limitations  in  the  numerical  solution 
of  the  radiation  transport  equations,  which  are  still  not  completely  resolved. 
Questions  as  to  the  treatment  of  thin  (optically)  hot  zones  adjacent  to  thick 
cold  zones,  shock  fronts,  radiative  fronts,  the  use  of  Planck  emission 
means  or,  for  that  matter,  what  is  the  appropriate  average,  etc. ,  still 
remain.  The  frequency-dependent  calculations  lead  to  questions  regarding 
the  selection  of  frequency  groups  and  the  use  of  various  frequency- group 
averages.  The  present  subroutines  allow  for  a  basic  solution  of  the  spherical 
transport  equations,  in  which  many  improvements  relative  to  the  above 
questions,  as  well  as  additions  in  respect  to  such  matters  as*  retardation, 
scattering,  etc. ,  can  be  made. 
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SECTION  H 


NUMERICAL  SOLUTION  OF  THE  TRANSPORT  EQUATION 


The  radiation  routines  described  herein  contain  a  formulation  based 
on  numerical  solution  of  the  radiation  transport  equation  along  a  selection 
of  sampling  rays  through  the  sphere.  Relevant  averages  over  the  angular 
distribution  are  obtained  by  numerical  quadrature,  as  described  in 
Section  2.  3,  and  the  numerical  solution  of  the  transport  equation  along  the 
photon  ray  is  presented  in  Section  2.  1.  Criteria  for  selecting  the  sampling 
rays  are  discussed  in  Section  2.  2.  All  of  the  derivations  of  this  section 
apply  to  photons  of  a  particular  frequency;  integration  over  frequency  is 
discussed  in  Section  IV. 

The  radiation  transport  equation  in  spherical  geometry  that  describes 
the  changes  in  the  specific  intensity  Iv  of  photons  of  frequency  v  resulting 
from  pure  absorption  and  emission  according  to  the  local  thermodynamic 
equilibrium  assumption  is 


where 


B 


2h 


2  hv/0  , 
c  e  -  1 


,  .  n  -We. 

<rv  -  %(1  '  e  > 


and  crv  is  the  pure  absorption  coefficient.  The  scattering  coefficient  is 
assumed  to  be  negligibly  small  compared  to  the  absorption  coefficient. 
Additionally,  the  retardation  of  the  photons  is  neglected,  as  is  valid  when 
the  radiation  energy  is  small  and  temperatures  change  slowly.  The  resulting 
equation  describes  the  quasi-steady  intensity  field  resulting  from  the 
distribution  of  sources  existing  at  a  particular  time.  This  equation  is 
simplified  by  introducing  two  new  independent  variables,  x  and  y,  in  place 
of  r  and  p,  where 


y  = 


x  =  r  n. 
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r 


(2.1) 


The  distance  Wong  fee  photon  ray  is  measured  by  x  and  fee  distance  of 
closest  approach  of  fee  undeviated  ray  is  given  by  y.  Since  each  ray  is 
characterized  by  a  particular  y-value,  it  is  also  called  a  y-lice.  These 
relations  are  illustrated  in  Fig.  2.  1,  where  fee  geometry  of  two  y-lines 
is  shown  in  relation  to  fee  radius  and  polar  angle  (Fig.  2.  la)  and  these 
quantities  are  translated  into  the  corresponding  two  rays  at  a  typical  point 
on  the  spherically  symmetric  surface  (Fig.  2.  lb).  In  terms  of  x  and  y, 
the  transport  equation  becomes 


~~  =  <r'  (B  -  I  ) 
8x  v  v  v 


an  equation  which  could  equally  well  be  considered  the  fundamental  defining 
equation  for  fee  intensity  in  that  it  describes  the  processes  that  change  the 
intensity  along  the  ray.  Defining  the  monochromatic  optical  depth,  t,  as 


(2.3) 


the  differential  equation  can  then  be  integrated  between  two  points,  t.  ^  and 
Tj,  on  a  ray: 


This  equation  forms  the  starting  point  for  the  numerical  solution  of 
the  transport  equation.  For  a  particular  ray  having  a  y-value  chosen 
according  to  the  prescriptions  given  in  Section,  2.2,  the  intensity  is 
evaluated  at  selected  points  along  the  ray,  starting  from  the  outside  boundary 
of  the  sphere  with  a  prescribed  boundary  value.  The  immediate  problem, 
then,  is  to  prescribe  the  method  for  approximating  the  quadrature  in  Eq.  (2.4). 
This  prescription  and  the  task  of  prescribing  boundary  conditions  are 
discussed  next. 


2.  1.  INTEGRATION  ALONG  Y-LINES 

In  the  finite-zone  calculation  of  the  SPUTTER  code,  the  value  of  the 
source  function,  Bv,  is  known  only  at  a  series  of  radii  corresponding  to 
the  average  positions  of  zones  in  the  calculation.  Intensities  are  needed, 
however,  at  all  interfaces  between  zones  in  order  to  determine  radiation 
fluxes.  Consequently,  it  13  necessary  to  construct  an  interpolation  function 
between  known  values  so  that  the  integration  of  Eq.  (2.  4)  can  be  carried  out. 
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Of  the  many  alternatives,  a  simple  scheme  is  chosen  which  permits  the 
integration  to  be  performed  in  closed  form. 

2. 1. 1,  General  Formula 

As  indicated  in  Fig.  2. 2,  the  source  function  is  usually  assumed  to 
be  linear  in  r  between  given  values  of  at  the  midpoints  of  zones 
t-+x.  At  certain  interfaces,  however,  when  die  criteria  described  in 
Section  V  are  satisfied,  the  source  function  is  assumed  to  be  constant  at 
fee  zone  value,  which  is  also  shown  in  Fig.  2.2.  Values  of  the  source 
function  B-  are  first  obtained  by  interpolation  at  fee  zone  interfaces  t. 


Fig.  2.  2 — Illustration  of  discontinuous  source  function 


■which,  together  wife  fee  midpoint  values,  form  fee  termini  for  the  integral 
of  Eq.  (2,4).  More  specifically,  in  fee  interval  t-_  j  <  t  <  t-,  the  source 
function  takes  fee  values 


where 


and 


B  =  a  +  b  r,  t .  <  t  <  T. 

-  1-1  —  —  X 


a_  = 


T 


i-1 


B  =  a+  +  bfr, 


T.  !  <  T  <  7.  , 

l-i~  -  1 


(2,5) 
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where 


B.  -  B.  i 


For  the  case  of  a  constant  or  step-fraction  source,  the  source  fraction  B 
takes  a  value  dependent  on  which  interface  of  the  zone  is  affected.  If  die 
left  interface  (t  =  t-_  j)  satisfies  the  criteria  for  a  constant  source, 

B  =  B.  j,  r  <  t  <  t.  ,  , 

If  die  right  interface  (t  =  r?)  satisfies  the  criteria, 

B  =  B.  ,,  t.  ,  <t<t.  . 

1-i.  1— ~  1 

1  2  4  2  * 

The  integral  of  Eq.  (2. 4)  can  be  evaluated  with  die  interpolation  function 
of  Eq.  (2.  5)  to  give  for  the  intensity 


I. 

i 


i-1 


+  Je 


-A/2 


1-2 


(2.6) 


where 


or.  i  =  a  +  b  (t.  -  1)  , 

i-j  +  +  a 

pi-i  =  a--a+Mb--V^-f  -  ‘)  ■ 

y.  i  ~  b  (1  +  A  -  t.)  -  a 
i-T  -  i 


In  these  expressions,  A  =  t j  -  The  coefficients  of  Eq.  (2.6)  can  be 

re-expressed  by  using  the  definitions  of  Eq.  (2.  5): 


a.  i  =  B.  - 


1-7 


P,.l 

*  2 


B.  -  B.  1 
i  1-2 

A/2 


B.  -  B.  i  B.  i  -  3.  , 

1  l“2  1-2  _ I"  1 


A/2 


A/Z 
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The  terms  ic  Eq.  (2.  7}  may  be  interpreted  as  containing  combinations  of 
numerical  approximations  to  die  values  of  the  source  function  and  the  r 
derivative  of  the  source  function  at  the  boundaries  of  die  interval. 


This  form  of  the  equation,  in  fact,  can  be  obtained  in  another  way 
starting  from  Eq.  (2. 4j.  Two  successive  integrations  by  parts  transform  die 
expression  for  Xj  into  die  following  equivalent  form: 


■.-(-SH- •  <?•!?), ,W 


a  b  -(t.  -r) 

— -  e  i  dr 


(2.8) 


in  terms  of  values  of  die  source  function  and  the  first  two  derivatives  of  the 
source  function  with  respect  to  r. 

Id  an  optically  thin  interval,  die  most  important  contribution  arises 
from  the  terms  I-  j  and  B,  which  represent  the  transmitted  intensity  and 
the  emission  from  the  zone.  The  derivative  terms  cancel  in  this  approx¬ 
imation;  this  is  perhaps  more  directly  indicated  by  Eq.  (2. 4),  In  die 
optically  thick  interval,  which  is  die  extreme  opposite,  only  the  first  two 
terms  evaluated  at  i  are  usually  of  significance.  The  terms  from  i  -  1  are 
strongly  attenuated  and  32B/3t2  in  die  integral  is  usually  small.  In  the 
limit,  the  diffusion  approximation  results  from  the  term  SB/3t)..  Between 
limits,  it  is  necessary  to  consider  the  integral  term  in  Eq.  (2.  8). 

If  A  is  not  too  large,  a  representative  mean  value  of  the  exponential 
in  the  interval  may  be  taken  to  give  for  the  integral  of  Eq.  (2.  8) 

f  £f«-<VT)— -A/2[n).-f).l  ■ 

"t.  ,  3t  L  '  1  '  l- 1-> 


and  thus  die  expression  for  intensity  becomes 

(» -£}■ *  [»).  -“) ,.,])••“  «.» 

This  expression  has  just  the  form  of  Eqs.  (2.  6)  and  (2,  7)  when  the  difference 
expressions  are  identified  with  the  derivatives. 

It  is  clear  from  the  derivation  of  Eq.  (2.  6)  that  the  resulting  intensity 
is  a  positive  quantity.  With  positive  values  for  zone  source  functions,  the 
linear  interpolation  expression  assures  that  the  integral  contribution  is 
always  positive.  Since  the  boundary  intensity  is  always  a  positive  quantity, 
the  positivity  of  all  intensities  is  assured. 


In  the  diffusion  approximation  limit*  however,  Eq.  (2.  9)  is  to  be 
preferred  over  Eqs.  (2. 6)  and  (2. 7),  since  in  this  limit  only  quantities  at 
interface  i  will  survive*  and 


I.  =  B. 
x  1 
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which  can  be  evaluated  more  accurately  than  the  corresponding  difference 
approximation  to  the  derivative  of  Eq.  (2.7).  The  point  is  that  B  depends 
only  on  radial  position  and  not  on  angle*  sc  that  an  interpolation  formula 
using  t  is  more  artificial  than  one  based  on  a  radial  quantity*  such  as 


3B  3B 
3t  "  M  dh  ' 


(2.10) 


where 


The  independent  variable  h  depends  only  on  r,  so  that  angular  integrations 
of  I  j  can  be  performed  explicitly  in  the  diffusion  approximation,  which  takes 
account  of  the  dependence  on  angle  of  Eq.  (2.  10).  A  difference  approximation 
can  also  be  based  on  this  expression,  assuming  that  B  is  linear  in  h,  i.  e. , 


(2.11) 


where  fi.  =  X|/r;  is  the  cosine  of  the  angle  of  the  ray  at  the  interface  with 
radius  r..  The  corresponding  equation  for  the  intensity  is  Eq.  (2.  6), in  which 


a.  i  ~  B.  -  p. 


B.  -  B.  i 
1  1-2 
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In  spherical  geometry,  however,  p  varies  along  the  ray.  so  the 
difference  approximation  in  Eq.  (2. 11)  is  not  identical  with  the  linear-in-T 
difference  approximation.  Consequently,  Eq.  (2. 12)  is  different  from 
Eq.  (2.  7).  Although  Eq.  (2. 12)  is  superior  in  the  diffusion  limit  in  that 
fee  surviving  terms  give  fee  correct  diffusion  expression,  unfortunately, 
for  fee  general  case,  fee  positivity  which  applies  to  Eq.  (2. 7)  is  lost; 
consequently,  it  is  not  clear  which  expression  is  best.  {In  slab  geometry, 
where  fi  does  not  change  along  fee  ray,  the  difference  equation  analogue  of 
Eq,  (2. 1C)  is  exact  and  the  positivity  of  the  equation  resulting  from  the 
linear-in-h  assumption  is  assured. )  The  greatest  danger  of  negative 
intensity  comes  from  regions  of  small  optical  depth;  however,  these  are 
calculated  by  a  limiting  form  of  the  equations,  as  described  below. 

Equations  (2,6)  and  {2. 12)  are  solved  in  fee  STRANS  subroutine  further 
described  in  Section  V.  Should  negative  intensities  result  from  the 
calculation,  they  are  replaced  by  zero. 

2.1.2.  Special  Calculation  near  x  =  0 

The  source  function  profile  shown  in  Fig.  2.  2  fails  to  take  inro  account 
all  of  the  information  available  in  fee  neighborhood  of  x  =  0.  Since  fee  source 
function  along  the  y-line  displays  symmetry  about  x  =  0,  it  is  clear  that 
its  derivative  should  be  zero  there,  a  condition  violated  by  the  preceding 
interpolation  rule.  Consequently,  a  special  calculation  is  performed  in 
fee  two  intervals  adjoining  x  =  0,  which  are  formed  by  the  double  penetration 
of  fee  zone  by  the  y-line,  as  illustrated  in  Fig.  2.3a.  In  the  half- intervals 
immediately  adjacent  to  x  =  0  {see  Fig.  2.  3b),  the  source  is  taken  to  be 
quadratic  in  x,  the  constants  of  which  arj  determined  by  requiring  that 
fee  derivative  be  zero  at  x  =  0  and  that  the  function  take  on  the  known  source 
value  at  the  center  of  the  zone.  The  source  in  the  remaining  intervals  is 
formed  from  fee  linear-in-x  assumption.  Since  the  interpolation  function 
is  contained  within  a  single  zone,  the  function  in  Fig.  2.  3b  has  exactly  the 
same  form  as  a  function  of  r  as  it  has  as  a  function  of  x. 

The  quantity  xi  is  the  x-coordinate  of  the  center  of  the  zone  having 
inner  radius  Tq  and  radial  thickness  Ar  (see  Fig.  2.3a), 

xj  =  [Ar{r  +  0.  25  Ar)]  2 
2  u 

and  corresponds  to  the  value  of  the  source  function  Bi.  Values  of 
the  source  corresponding  to  x  =  0  and  x  -  x^  are  and  Bj. 


iO 


Fig.  2.  3- -Assumed  spatial  dependence  of  source 
function  near  x  =  0 


To  illustrate  the  derivation  of  the  intensities  Iq  =  I(x  =  0)  and 
Ij  =  I(x  =  xj),  the  following  steps  in  the  calculation  of  the  source  integral 
for  the  interval  -Xj  <  x  <  0  are  indicated.  In  this  interval  the  source 
function  is  represented  as  two  parts  (see  Fig.  2.3b) 


B  -  a  +  b  x  ,  -x  <  x  <  -xj.  , 


where 


a  = 


xiBI  *  xIBi 
1  _2 _ 2 _ 1 

‘I 


and 


B i  ~  B, 


b  = 


B  =  a,  +  b, x  , 
+  r 


x  -  x, 
1  2 


•x,  <  x  <  0  , 
2 


where 


In  terms  of  these  quantities,  the  integral  of  Eq.  (2.4)  is  then 
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b 

+  —  [(A+l)e 
<r 
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“A  i_  *Ai 

-  (Aj+l)e  3]  +  a  (I  -  e  2) 

2  + 

b+  2  “Ai 

+  -r-  [2  -  (Ai  +  2At  +  2)e  2]  , 

n~  2  2 


where  Ai  =  <rx±  and  A  =  mcj.  Substitution  and  reduction  gives  the  following 
expression  for  1^: 


-A, 


-A 


X0  =  B0  +  F  +  6  2  [G  “  (Ai  +  1)F j  +  e  “  [K-Xj)  -  Bj  -  G]  , 


where 


2(Bi  -  B  ) 

F  - - 2 - - 

a  2 
Ai 
2 


and  G  = 


B,  -  Bi 
1  2 

A  -  Ai  ‘ 
2 


(2.13) 


In  an  analogous  fashion,  Ij  is  obtained  by  integration  over  the 
interval  0  <  x  <  x ^ , 


-(A  -  AO 


I  =  B  -  G  +  e 


-A 


[G  -  (1  -  Ai)F]  +  e  [IQ  -  BQ  -  F] .  (2.14) 


Equations  (2.  13)  and  (2.  14)  are  evaluated  in  the  "top  slice”  portions  of  the 
code  when  the  y-line  integration  reaches  x  =  0. 


2.1.3.  Small-optical-depth  Expansion 

If  the  optical  depth  is  very  small,  the  intensity  expression  in  Eq.  (2.4) 
takes  a  much  simpler  form, 


I. 
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7  B.  +  7  B. 

4  i  4  l-l 


+K* 


(2.15) 


{ 


r 


( 


12 


Kim. 


Although  this  result  is  the  limiting  form  of  Eqs.  (2.  6)  and  (2.  7),  but  not  of 
Eqs.  (2.6)  and  (2.  12),  the  terms  must  cancel  through  second  order  in  an 
expansion  in  A  before  the  first  surviving  term,  derived  in  part  from  the 
quadratic  terms  of  the  exponentials,  is  obtained.  Consequently,  for 
sufficiently  small  argument,  the  finite  number  of  figures  used  in  the 
exponential  will  render  the  result  inaccurate.  For  the  exponential  from 
the  IBM-7044  system,  this  restricts  the  argument  to  a  number  greater 
than  ~2x  10“^;  but  with  the  lower-accuracy  fast  exponential  (see  Section  V), 
the  argument  must  be  somewhat  larger.  Since  the  relative  error  approx¬ 
imately  equals  the  argument  of  the  exponential,  the  criterion  for  using 
Eq.  (2.  15)  in  the  STRANS  subroutine  is  now  set  at  A  <  2  x  10“^.  With 
this  value,  the  greatest  relative  error  arising  from  the  expansion  and 
cancellation  should  be  on  the  order  of  1  percent. 


For  the  quadratically  interpolated  intervals  near  x  =  0,  the  small- 
optical-depth  result  is  different  from  Eq.  (2.  15)  as  follows: 


■  "-v +  [i B. +  7  V  +  (I  Bo  - 1 B1  -  i  BiK 


(2.16) 


The  same  criterion  for  performing  the  calculation  of  Eq.  (2.  16)  rather 
than  that  of  Eqs.  (2.  13)  and  (2,  14)  is  used,  i.  e.  ,  A  <  2  x  10“^. 


2.1.4.  Y  =  0  Ray 

The  special  case  of  the  ray  passing  through  the  center  of  the  sphere 
and  having  y  =  0  is  required  (as  described  below)  in  each  use  of  the  STRANS 
subroutine.  Equation  (2.  12)  is  still  appropriate  for  this  case  but  in  a 
simplified  form.  For  inwardly  directed  photons,  n  =  -1;  for  outwardly 
directed  photons,  n  -  +1.  A  separate  section  of  STRANS  is  used  for  the 
y  =  0  calculation  in  order  to  simplify  the  code  and  to  take  into  account  that 
the  angular  integration  need  not  be  performed. 

2.1.5.  Boundary  Conditions 

Integration  of  the  transport  equation  to  obtain  intensities  is  performed 
through  the  thickness  of  a  spherical  shell,  called  a  "trans"  region.  At 
intersections  of  y-lines  with  the  inner  and  outer  surfaces  of  each  shell  it 
is  necessary  to  supply  the  starting  value  of  the  intensity  I^_j  required  in 
Eq.  (2.6).  Three  classes  of  boundary  conditions  occur: 


1.  The  trans  region  outside  boundary  coincides  with  outside 
radius  Rjg+j  of  the  SPUTTER  calculation  and  a  prescribed 
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function,  Iq,  is  applied  as  an  outer  boundary  value: 


I<Rm+i*  ^  =  V  M  -  0  • 


<2.17) 


The  trans  region  inside  boundary  is  at  the  center  of  the 
sphere,  R  =  0. 

2.  The  center  of  the  sphere  provides  a  boundary  condition  for 

the  y  =  0  ray  at  R  =  0  in  the  trans  region.  At  R  =  0,  the  intensity 
is  isotropic,  giving  the  condition  for  starting  the  outward-directed 
ray  calculation: 


I(R  =  0,  M  =  +i)  =  I(R  =  0,  ft  =  -1)  .  (2.18) 


3.  All  other  trans  boundaries  are  bounded  by  regions  in  which  the 
diffusion  approximation  is  valid  (see  Section  III).  Consequently, 
the  boundary  surface  intensities  on  contiguous  trans  regions 
inside  or  outside  of  a  diffusion  region  are  given  by  the  diffusion 
approximation  intensity  derived  in  Section  III; 


8B\ 

Mi-1  ah/i-1 


(2.19) 


2.  2.  Y-LINE  POSITIONING  REQUIREMENTS 


The  correct  placement  of  y-lines  through  regions  which  are  changing 
rapidly  in  temperature  and/or  optical  depths  is  important  for  an  economical 
and  accurate  solution  of  the  transport  equations  using  the  direct  integration 
method.  The  advantage  of  being  able  to  locate  y-lines  only  where  they  are 
needed  is  a  decrease  in  computational  time  by  reducing  the  number  of  y-lines 
calculated. 

The  decision  as  to  where  to  place  a  group  of  rays  depends  on  the 
radial  position  of  the  region  in  question.  The  inner  core  (or  zone)  will  be 
the  yjast  resolved,  but,  in  general,  this  io  a  region  of  isothermality. 

For  jjxample,  in  fireball  calculations,  there  is  a  central  diffusion  region 
bounded  by  a  region  having  rapidly  changing  conditions  due  to  radiative 
fronts  and  hydrodynamic  shocks.  Under  these  conditions,  there  are  fewer 
y-lines  placed  through  the  central  core,  many  through  the  transition  regions, 
and  none  outside  the  zones  containing  sources. 

The  following  description  of  our  placement  criteria  in  SRADTN  and 
STRANS  is  based  on  these  considerations. 
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Z.  L.  i.  X-storage  Limitation  and  Frequency-independent  Criteria  in 
SEADTN 


Jxc  order  to  remove  toe  repetitious  calculation  ox  x-values  at  inter¬ 
sections  of  the  y-liaes  with  radii  x  -  * Jr ^  -  Y^,  a  storage  array  is  set  up 
for  fee  x*s  as  weil  as  for  fee  quantity  y2  (used  in  STRANS).  The  size  of 
this  array  is  limited  by  available  storage;  fee  present  limit  of  24C0  in  the 
x-block  allows  for  approximately  70  y-lines.  The  principal  criterion  for 
selecting  a  set  of  y-lines  to  be  used  in  STRANS  is  to  place  y-lines  at 
boundaries  where  fee  temperature  is  changing.  For  example,  if  log  (9j+j /8^) 
is  greater  than  CVB  (CVB  =  0  gives  the  maximum  number  of  y-lines),  a 
y-line  is  placed  at  feat  boundary,  but  if  fee  gradient  of  log  0  is  less  than 
CVB,  feat  zone  is  skipped.  The  practical  use  of  this  criterion  will  result 
in  reducing  fee  calculational  time  for  y-lines  inside  isothermal  regions. 

No  y-lines  are  used  In  regions  outside  zones  with  fee  temperature  less  than 
0.  05  ev  (0. 025  ev  is  the  temperature  for  ambient  air). 

2.2.2.  Transport  Criteria  per  Frequency  Group  in  STRANS 

In  fee  transport  subroutine  (STRANS),  y-lines  from  fee  set  established 
in  SRADTN  are  selected  using  criteria  based  on  the  frequency-dependent- 
source  and  optical-depth  gradients.  If  a  diffusion  region  having  an  outer 
boundary  rp  exists  inside  fee  transport  region,  y-lines  are  placed  as  near 
as  possible  to  0.  5  rp,  0.  75  r^,  and  rp  penetrating  fee  diffusion  region.  In 
addition,  if  inside  of  fee  transport  region  the  gradient  inequalities 


(b  e4)  i  -  (b  e4)  i 

_J _ *+*  J _ lZ± 

(bi04)i+i  +  <V4)i  ± 
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are  satisfied,  where  the  input  number  GL  is  usually  around  0.3,  a  y-line 
is  added  at  the  interface  i.  Finally,  in  isothermal  source  regions  outside 
the  diffusion  core,  every  fifth  ray  is  selected.  These  selection  rules  could 
be  refined,  but  at  the  present  time  they  afford  a  reasonable  representation  lo 
be  used  in  the  angular  integration  of  the  intensities. 
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2. 3.  ANGULAR  INTEGRATIONS 


Intensities  calculated  in  the  integrations  along  y- lines  are  used  in 
several  ways  in  the  SPUTTER  calculation*  In  addition  to  displaying  die 
intensities  in  a  special  edit  (see  Section  5.  5),  integrals  over  angle  of  die 
intensity  are  used  in  the  calculation  cf  the  evolution  of  the  system  in  time. 
Additional  integral  quantities  are  of  interest  in  evaluating  the  progress  of 
the  problem  and  in  preparing  thermal  environment  data  for  subsequent 
calculations.  The  desired  integral  quantities  are  formulated  in  Section  2.  3.1. 

To  calculate  die  desired  integral  quantities  from  the  intensities 
obtained  in  the  integrations  along  a  discrete  sampling  of  y~  lines,  as 
described  in  Section  2.  1,  it  is  necessary  to  perform  numerical  quadratures. 
At  a  given  interface  radius  no  regular  interval  in  the  angle  variable 
(either  fi,  x,  or  y)  is  available.  Consequently,  die  numerical  quadrature 
is  based  on  an  interpolation  expression,  described  in  Section  2.  3.  2,  in 
which  no  regularity  can  be  assumed.  Finally-  in  Section  2. 3. 3,  the 
numerical  quadrature  formulas  are  given. 

2.3.1.  Integral  Formulations  of  flux.  Energy,  and  Pressure 

One  of  the  most  important  of  the  angular  integrals  of  intensity  is  the 
radiative  flux,  the  divergence  of  which  gives  the  radiative  contribution  to 
the  rate  of  ckange  of  the  material  energy.  The  net  flux  is  a  vector  quantity, 
the  component  of  which  in  the  direction  of  the  unit  vector  r  j,  4>i)  (i» 
ergs/cm^-sec),  is  given  by  the  expression 


4»{^1>  =  c  J  I  cos  *  da  ,  (2.2G) 

which  also  depends  on  the  spatial  position  r. 

In  Eq.  (2.  20)  the  intensity  is  a  function  of  position  r  and  direction, 
denoted  by  the  unit  vector  r,  the  Cartesian  components  of  which  can  be 
expressed  in  terms  of  the  polar  angle  9  between  r  and  the  radial  direction 
k  and  the  azimuthal  angle  4>  measured  between  the  meridian  planes  of  r 
and  an  arbitrary  direction  j  normal  to  It,  and  hence 


r  =  cos  0  k  +  sin  0  (sin  4*  +  cos  <p  j)  • 


The  angle  between  r  and  is  ®,  so  that  cos  ®  =  r  •  r^  , 
of  solid  angle  is  dfl  =  d$  sin  0  d0.  The  net  flux  in  the  r’j 


and  the  element 
direction  is  given 
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by  the  integration  of  d2  over  fee  entire  sphere.  The  forward  arsd  backward 
currents,  and  4",  however,  are  obtained  by  integration  over  fee  hemi¬ 
spheres  H*  and  H“  about  and  about  -r,  ,  respectively;  i.  e.  , 


I  cos  •  dQ 


9 


9  (r,) 

a. 


I  cos  •  d2  . 


i2.2i) 


In  terms  of  these  quantities,  fee  net  flux  is 


♦(r^)  =  {r j)  -  $  (Tj)  .  (2. 221- 

In  spherical  geometry,  in  which  I  depends  on  6  but  not  on  6,  fee 
integration  over  azimuthal  angle  can  be  performed,  and  fee  result  depends 
only  on  fee  radial  distance  r  and  fee  component  direction  rj .  hi  order  to 
evaluate  fee  progress  of  fee  calculation  and  to  provide  relevant  thermal 
environment  data  for  subsequent  calculations,  three  representative  currents 
are  formed- -fee  forward  and  backward  currents  in  fee  radial  direction  and 
fee  current  in  fee  direction  normal  to  fee  radius  vector.  Only  one  such 
lateral  current  is  needed  since  the  currents  in  all  lateral  directions  are  fee 
3ame  for  spheres.  For  the  forward  current,  <?T,  fee  quantities  in  Eq.  (2.21) 
are  rj  =  k,  cos  •  =  cos  0  =  fi,  and  fee  integral  is  over  0  <  9  <  x/2  and 
0  <  4  <  2»;  feus 


(2-23) 


For  the  backward  current,  4*  »  fee  quantities  are  r^  -  k,  cos  €>  =  n,  and 
the  integral  is  over  w/2  <  9  <  ir,  0  <  $  <  2x;  thus 


4> 


-2trc 


Ip  d n 


(2.24) 
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For  fee  lateral  current,  t  r  the  forward-current  expression  of  £q.  (2.  21) 
may  be  evaluated  by  taking  r-  =  i,  cos  •  =  JT-jJ  sin  and  0  <  0  <  », 
so  shat 


.0 

♦ 


(2.25) 


la  spherical  geometry  fee  term  entering  fee  energy  equation,  V  -  9, 
in  which  $  i>  fee  radiative  net  flax  vector,  takes  fee  form 


i_  JL 

2  3r 


in  which  is  fee  net  flux  component  in  fee  radial  direction  k.  Two  of  fee 
above  quantities  are  related  to  9r  in  feat 


(2.  26) 


Consequently,  of  fee  three  quantities,  9r»  9  ,  and  9  ,  only  two,  9r  and  6  , 
are  formed. 

Two  additional  angular  integrals  cf  the  intensity  are  useful  in  evaluating 
the  progress  of  the  transport  calculation  and,  in  addition,  they  are  used  In 
more  accurate  formulations  of  the  radiation  transport  equation.  The  first 
of  these  quantities,  the  radiation  energy  density,  Er  (in  ergs/cm^),  gives  a 
quantitative  measure  of  the  energy  stored  in  the  radiation  field  for  com¬ 
parison  with  the  material  energy.  If  Er  is  not  negligible  (as  assumed  to  be 
the  case  in  this  formulation),  it  should  be  taken  into  account  in  the  radiation 
transport  equation  in  which  retardation  is  included  (see  Section  2.  1  of 
Volume  VI).  The  same  angular  integral  also  plays  a  role  in  the  Thomson 
scattering  integral  of  the  Compton  scattering  by  free  electrons,  as  discussed 
in  Section  3.  1  of  Volume  VII  (the  effect  of  which  is  also  neglected  in  this 
formulation). 

The  radiation  energy  is  given  by 


which,  for  spherical  geometry,  in  which  the  azimuthal  integration  can  be 
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performed,  redaces  to 

(2-27) 

» 

The  second  quantity  is  the  radiation  pressure,  Pg(rj)  {in  ergs/cm  ), 
which  is  defined  as  the  net  rate  of  transfer  of  the  radiant  momentum  com¬ 
ponent  in  the  r  j  direction  across  the  unit  surface  whose  normal  is  also  in 
the  rj  direction.  In  terms  of  cos  •  =  r  j  •  r,  fee  cosine  of  the  angle  which 
the  photon  beam  makes  wife  fee  surface  normal,  fee  pressure  exerted 
across  rj  by  fee  radiation  field  is 

pR(ri>=  f  Icos2«dO  , 


in  which  fee  integral  over  solid  angle  extends  over  fee  complete  sphere. 

The  pressure  Pjj{rj>  is  also  obtained  from  fee  radiation  pressure 
tensor  P  by  the  operations  r  j  -  P  •  r, .  For  spherical  geometry  in  which 
fee  azimuthal  angle  integration  can  be  performed,  all  of  fee  off-diagonal 
elements  of  fee  pressure  tensor  vanish,  whereas  all  of  fee  diagonal  elements 
can  be  expressed  in  terms  of  fee  diagonal  element  in  fee  radial  direction, 

1 

PR  (k)  =  2 v  J  I  /  d«  s  Pr  ,  ( 2 .  28) 

and  fee  radiation  energy  integral  of  Eq.  (2.  27).  Furthermore,  the  integral 
of  Eq.  {2.  28)  is  also  contained  in  the  expression  for  fee  Thomson  limit  of 
fee  Compton  scattering  into  the  beam  in  both  plane  and  spherical  geometries. 

In  summary,  the  radiation  integrals  are  calculated  in  the  radiation 
subroutine  STRAN5  in  terms  of  the  following  quantities:  <j>+,  Eq.  (2.23); 

4>°,  Eq.  (2.25);  <f»r,  Eq.  (2.26);  ER,  Eq.  (2.27);  and  PR,  Eq.  (2.28). 

2,3.2.  Angular  Interpolation  of  Intensities 

To  form  numerical  approximations  to  fee  angular  integrals  derived 
above,  it  is  necessary  to  select  a  quadrature  formula.  Equivalently,  the 
rule  must  be  specified  for  the  angular  interpolation  of  intensities  at  a 
given  radius  between  values  calculated  at  different  values  of  x  (or  n)  in  the 
integrations  along  y-lines.  Unfortunately,  from  the  point  cf  view  of  the 
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angular  integration,,  as  a  result  of  the  application  of  the  y-line  selection 
criteria,  the  values  of  x  arc  sot  regularly  spaced  so  no  simple  high-order 
rule  is  easily  applicable.  *In  the  current  version,  the  interpolation  is 
performed  within  the  interval  between  adjacent  x  values  by  making  use  of 
the  two  values  of  intensity  at  the  end  points  of  the  interval  only. 

An  accurate  integration  will  result  if  the  interpolation  formula  conforms 
to  the  actual  angular  dependence  of  the  intensity.  To  illustrate  the  range  of 
possible  dependencies  which  may  occur  in  a  fireball  calculation.  Figs.  2. 4 
through  2. 7  snow  a  selection  oi  curves  of  angular  dependence  of  intensity 
corresponding  to  selected  radii  and  photon  frequencies.  It  is  clear  from 
these  curves  that  no  simple  dependence  will  be  universally  successful.  In 
fact,  these  curves  point  up  the  desirability,  in  a  mow,  advanced  code,  of 
having  an  interpolation  formula  in  which  information  from  several  angular 
positions  is  used. 

In  the  preparation  of  the  current  version  of  the  STRANS  subroutine, 
two  interpolation  rules  were  investigated:  linear-in-p2  and  linear -in-p 
dependence.  Although  the  ft2  interpolation  is  suitable  for  some  cases, 
such  as  that  of  Fig.  2. 6,  fee  linear-in-/*  dependence  approximates  fee 
behavior  of  most  of  fee  cases  illustrated.  Furthermore,  both  the  diffusion 
approximation  and  fee  transparent  isothermal  regime  have  this  dependence 
as  a  limiting  value.  Consequently,  fee  linear-in-p  interpolation  formula 
has  been  incorporated  into  the  STRANS  routine. 

2.3.3.  Numerical  Quadrature  Formulas 

The  contribution  of  fee  angular  interval  (pj,  Mjxj)  to  these  integrals  at 
position  r  is  derived  in  the  linear-in-p  approximation.  The  intensity  in  this 
approximation  is 


I  =  a  +  bp,  /*i<M<Mi+1  » 


where 


a  =• 


ft.,  .1.  -  ft.  I.,  .  x  I.  -  x.  I 

l+l  i  i  i+I  i+l  i  i  i+1 


**i+I  '  “i 


X.,  .  -  X. 
l+l  1 


b  ~  **  _  r_j±I_L± 
=  #1i+i  •  =  S+i  ■  xi 


I.  is  the  value  of  the  intensity  at  According  to  the  y-line  integration 
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scheme,  and  M^+i  correspond  to  the  angles  formed  at  the  intersections 
of  two  adjacent  y-lines  with  the  interface  radius,  at  which  point  the  integral! 
are  to  be  formed.  Two  intervals  are  symmetrically  located  about  n  =  0  by 
a  pair  of  y-lines,  as  shown  in  Fig.  2.8.  It  is  clear  from  the  figure  that 


Fig.  2.  8 — Angular  intervals  between  adjacent  y-lines 


these  two  intervals  can  be  treated  together  in  the  integrations  over  the  interval 
-1  <  n  <  1.  Denoting  intensities  with  positive  n  by  I+  and  those  with  negative 
n  with  I“,  the  integrals  cam  be  grouped  according  to  whether  the  integrand 
is  even  or  odd.  For  even  integrands. 


1  1  , 

£  i  m  dM  =  J  a  + 


i )  m  d  n  , 


(2.29) 


where  f(-  |m  |  )  =  f(  |m  |),  and  for  odd  integrands, 


1  1 
J  i  m  a (i  - 


I  )  m  dM  , 


(2.  30) 


where  f(-  |/i  |)  =  -f{  jji  j),  I(  |m  |)  5  I  (m)  and  I(-  \n  j)  =  I~(/i)  in  which  0  <  n  <  1. 

For  a  given  radius  r,  the  quantities  p,  x,  and  y  associated  with  the 
interval  (|^.J,  i^+i  |)  satisfy  the  following  inequalities: 


>  KI’  xi+i>xi’  yi>yi+i  • 
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where  x.  =  lu.  Jr  and  y.  =  vr2  -  x?  . 

x  1  i 1  J  i  i 

In  terms  of  the  above  quantities,  the  contributions  to  the  integrals  of 
the  interval  { jp.J,  |AH+1  J)  in  the  linear-in-p  approximation  are  obtained  as 
follows: 


Forward  Flux,  <}>  (Eq.  (2.23)): 


Mi+i 


T+  .  a+  .  2  2.  b+  .  3  3, 

Ipdp  =  T(M.+  1  -#*t)+T  (Mi+1  “  V 


+  2*i>  <  +  <2*i+,  +  %»  I+i+i] 
6r 


(2.31) 


Net  Flux,  4>r  (Eg.  (2.  26)): 


Mi+1 


/  <I+  *  O  M  *  =  ^  <1  -  >1'  +  1  '  ^ 

Jp. 

x 


fa.,  .  '  X?)  +  - 

‘  2  J  K*»i  *  2xi><li  -  V 

6r 


(2.  32) 


+  (2*.+  1  +  *;>«;„  -  I;+i>] 


Lateral  Flux,  <{>  (Eq.  (2.25)): 


J  (I+  +  i")  7l  -  M2  dp  =  (a  "t,  j^x. 


i+1  Vl  •  Xiyi 


if  .  -1  Xi+1  .  -1  Xi  \ 
+  r  (sin  — - sin  —  J 


t  <b+  +  b:>  ,3 .  3  , 

,3  iyi  M+r 

3r 


(2.  33) 
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Radiation  Energy,  Er  (Eg.  (2.27)): 


M... 

r  1+1  +  _  + 

I  (I  +  I  )  d#i  =  (a  + 

**  f*. 


(^r  »  1  *  ^*)  1  J.  _ 

=  _iii — L  d+  + 1.  + 1+  +  i  )  . 

2  r  i  l  x+1  l+I 


(2.  34) 


Radiation  Pressure,  PR(Eq.‘  (2.28)): 


,‘ltl  •  2  .  _  (a+  +  a"}  ,  3  _  .3,  ,  (b+  *  bj  .4  4 

4  VMi+l  i 1 


f  +  -  2  (a  +  a  )  ,  3  3,  J 

/  (I  +  I  )/x  dM  - - 3 - l^i+i  ~  V  ' 


(2.35) 


=  ^-3  fd{  +  \)(y  -  4x13)  +  <4i+  W(4xm  '  Y)1  ’ 


12r 


2  2  2 

where  y  =  (2r  -  y.  -  y.+1)(x.  +  x.+  1). 
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SECTION  in 

THE  DIFFUSION  APPROXIMATION 


The  radiation  transport  equation  in  the  limiting  case  cf  an  optically 
thick  medium  admits  of  the  diffusion  approximation  in  which  the  expression 
for  the  radiation  intensity  is  greatly  simplified;  only  the  local  properties 
affect  the  radiation  intensity  at  the  point  in  question.  An  expansion  of  the 
radiation  source  function  Bv  about  the  point  r  permits  the  intensity  ipfyt}  of 
the  radiation  field  in  the  direction  making  an  angle,  whose  cosine  is  ft,  with 
the  radial  direction  to  be  formed. 


3.  1.  DIFFERENTIAL  FORM  OF  THE  DIFFUSION  FLUX 


The  general  solution  of  the  transport  equation  forms  the  starting  point 
of  the  derivation.  The  integral  expression  for  the  intensity  applicable  to 
all  geometries  is 

T 

I{t)  =  f  B(T')e'(T'T  }  dr'  , 

J- oo 


where  T  =•  SqV  kv p  ds,  in  which  kv  is  the  monochromatic  absorption  coefficient 
(in  cm^/g)  at  frequency  v.  By  expanding  B(t')  in  series  about  the  point  t, 
i.  e. , 


B(t')  =  B(t)  +|f  (t! 


.  I  82B  .  .  ,2  , 

T)  +  2  2  T  '  T)  + 
2  T 


the  intensity  becomes 


I  = 


B  - 


9B  i  92B 


or 


I  =  B 


JL  M  +  (L.  JL  fj_  8bN\ 

Kp  dr  Kp  dr  \/cp  8r  / 


for  spherically  symmetric  geometry. 
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The  diffusion  approximation  results  from  retention  of  only  the  first 
two  terms,  so  that  the  diffusion  intensity  is 


T  =  K 

*  jq>  2r 


(3.1) 


and  the  monochromatic  diffusion  flux  or  and  radiation  energy  E  are 

*  *  K 

_  r*  r  .  4tc  1 

*c  x,4‘4‘  =  *”^ 


as 

cp  0r 


(3.2) 


er  =  2’ 


I  djx  =  4rB 


3.  2.  CRITERIA  FOR  THE  SELECTION  OF  DIFFUSION  REGIONS 


The  criteria  for  the  validity  of  the  diffusion  approximation  can  be 
obtained  by  examination  of  the  above  derivation — namely,  that  the  expansion 
of  the  source  function  be  justified  and  that  the  expansion  converge  rapidly 
so  that  the  neglect  of  all  but  the  leading  terms  is  valid.  If  the  source  function 
is  linear  in  t*  at  the  point  in  question  and  is  also  linear  for  a  distance  of 
the  order  of  one  mean  free  path  on  either  side  of  the  point, the  criteria  are 
satisfied.  These  criteria  are  difficult  to  quantify  since  they  ref?r  to  a 
finite  region  containing  the  point  in  question.  If  all  of  the  terms  (or  a  large 
number  of  them)  were  checked  for  rapid  convergence,  this  would  imply 
(making  a  smoothness  assumption)  that  the  diffusion  criterion  is  met.  It  is 
not  possible  with  finite  differences,  however,  to  form  the  higher-order  local 
derivative  approximations. 

In  the  SPUTTER  subroutine  SRADTN,  criteria  designed  to  give  an 
indication  of  both  the  local  and  nonlocal  behavior  have  been  employed.  First, 
at  the  zone  interface  at  which  the  intensity  and  flux  are  to  be  evaluated,  the 
inequality 

lffl<<B  <3-3> 

is  required.  In  this  expression  h  =  St: p  dr  is  the  radial  optical  depth;  the 
derivative  is  approximated  by  the  centered  first  difference  of  B  between 
adjacent  zones.  The  resulting  expression,  of  course,  contains  some  non¬ 
local  aspects  resulting  from  the  finite  difference  approximation,  which 
ensures  that  when  neighboring  zones  are  optically  thick,  no  nonlocal  source 
perturbation  is  close  enough  to  invalidate  the  diffusion  approximation. 


However,  to  provide  for  the  cases  when  a  source  perturbation  is  located 
a  fraction  of  as  optical  depth  from  as  interface  meeting  the  condition  of 
Eq.  (3. 3),  the  diffusion  region  is  constricted.  Starting  from  the  closest 
interfaces  outside  the  diffusion  region  (where  Eq-  (3.  3}  is  sot  satisfied), 
all  of  those  interfaces  lying  within  a  prescribed  number  of  mean  free  paths 
are  removed  from  the  diffusion  region. 

The  criteria  used  in  SPOTTER  are  controlled  by  input  numbers.  The 
criterion  of  Eq.  (3-  3)  uses  the  input  number  HCB: 


|TGf  <  HCB  x  Y2 


(3-4) 


where  TG  is  the  difference  approximation  to  the  gradient  and  Y2  is  the 
source  function  evaluated  at  the  interface  by  interpolation.  The  second 
criterion  uses  the  input  number  HVB  (in  mean  free  paths).  If 


lQ3(I)  -  Q3(J)j>  HVB 


(3.5) 


then  the  interface  with  index  J  which  satisfies  Eq.  (3.  4)  is  removed  from 
the  diffusion  region.  In  Eq.  {3.  5),  Q3  is  the  radial  optical  depth  and  I  is 
the  index  of  the  nondiffusion  interface  adjoining  the  diffusion  region. 

Although  the  diffusion  calculation  is  considerably  faster  than  the 
transport,  the  establishment  of  two  transport  regions  separated  by  the 
single  zone  requires  still  more  calculation  to  set  up  y-lines  and  perform 
bookkeeping  operations.  To  avoid  the  duplicate  setup  calculations  required 
for  an  additional  transport  region,  a  test  is  made  to  eliminate  a  diffusion 
region  consisting  of  a  single  zone. 


3.  3.  DIFFERENCE  FORM  OF  THE  DIFFUSION  F  ^UX 


The  diffusion  intensity  derived  above  is 


I  =  B 


£_  SB 
icp  dr 


In  the  group  frequency  approximation  of  SPUTTER,  the  intensity  integrated 
over  a  frequency  interval  (v.,  v.  ,)  is  reauired: 

j  j+1 

Vj+1  3B  dv 


rvj'H  J  r 

j+I  „  ,  ^ 

804  r 

f  1  dv  =  J 

J  V  .  v 

B  dv  -  - 
p 

,  (CO 

9S4'v 
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In  terms  of  the  partial  Rcsseland  mean  absorption  coefficient 


/ 


8B 

Vj  394 

K.  1 - - - 

J  v- 


dv 


/ 


3B  dv 


_  4  x. 
v.  38  v 
J 


the  frequency  group  intensity  becomes 


v.  ,  » 

v/'— / 

■*  * »  it 


f 


Vi 


3B  ^ 

V.  ,  v.  ,  .  J  — 7  dv  . 

-  J+1  f  j+1  4  v  4 

-  (3.6) 

P  or  x. 

v.  v.  r  3 

j  ; 


It  is  desired  to  evaluate  this  quantity  at  each  zonal  interface  in  the 
mesh.  Since  the  known  quantities  are  the  zone  temperatures  and  densities, 
the  absorption  coefficients  Kj  and  the  integrated  source  functions  X6  =Jb  dv 
are  first  evaluated  not  at  the  interfaces  but  at  positions  representative  of 
each  zone. 

The  question  remains  as  to  how  best  to  approximate  the  derivatives 
and  interpolate  for  the  coefficients  in  Eq.  (3.  6)  at  the  interfaces  from  the 
quantities  available  at  zone  positions.  The  answer  depends  on  the  tem¬ 
perature  and  density  profile  across  the  interface  from  which  these  terms 
could  be  calculated  directly.  Since  the  profile  is  not  known,  we  must  select 
a  reasonable  approximation  which  will  permit  the  calculation  to  be  carried 
out.  In  fact,  the  appropriate  profile  depends  on  the  events  which  have  taken 
place  in  the  calculation  and  on  the  energy  transport  mechanisms  of  greatest 
importance  in  it.  As  extreme  examples,  a  problem  dominated  by  hydro- 
dynamic  might  have  quantities  determined  by  passage  of  a  strong  shock 
and  subsequent  linearization  in  mass  coordinates  of  the  pressure  behind 
the  shock,  whereas  a  radiation-dominated  diffusion  problem  is  characterized 
by  linearity  of  the  radiation  potential,  which,  in  turn,  depends  on  the 
Rosseland  opacity.  Of  course,  such  detailed  information  about  the  progress 
of  a  problem  is  generally  unavailable,  so,  at  best,  an  approximation  based 
on  over-all  accuracy  is  needed. 


Since  the  terms  under  consideration  are  the  radiation  diffusion  equa¬ 
tions,  the  interpolation  is  performed  in  a  way  to  give  greatest  accuracy 
when  the  diffusion  terms  are  most  important- -namely,  when  the  profile 
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is  being  determined  entirely  by  radiation  diffusion.  It  is  also  desirable  to 
redace  the  number  of  coefficients  requiring  interpolation.  This  can  be  done 
by  noting  the  identity 


3  f^1  _  .  364  r")*1  &B  , 

or  J  3r  I  _  4 

-'v.  Jv-  30 


and  by  forming  the  variable  r  =  S  P  Sj  dr.  In  terms  of  these  quantities,  the 
intensity  can  be  written  as 


■/ 


j+1 


•/ 


Vi 


B  dv 


B  dv  -ji 


3t 
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SECTION  IV 


FREQUENCY  INTEGRATION 


Equations  derived  in  Sections  II  and  III  which  are  applicable  to  a 
particular  frequency  of  the  radiation  field  are  of  limited  usefulness  in  the 
SPUTTER  calculations.  Although  in  principle  a  calculation  at  a  particular 
frequency  might  be  valuable  for  comparison  with  high -re solution  spectro¬ 
scopy,  in  practice  no  such  data  have  been  available.  Of  much  more  use  are 
intensities  averaged  over  a  wide  frequency  band.  These  quantities  can  be 
compared  with  data  from  wide -band  measurements  and,  most  important  of 
all,  can  be  summed  for  use  in  the  energy  integration  in  the  SPUTTER  code. 
The  quantities  to  be  summed  are  the  frequency-integrated  radial  flux 
component,  the  radiation  energy  density,  and  the  radiation  pressure.  For 
performing  interaction  calculations,  it  is  also  valuable  to  form  other 
components  of  the  radiation  flux. 

Basically,  the  quantity  which  is  required  for  each  of  the  above 
applications  is  the  frequency-group  intensity  I^j, 

v.  . 

fJ+1 

I..  =  I.  dv  .  {4.  1) 

lJ  K.  1 

J 

Then,  for  example,  this  quantity  can  be  integrated  over  angles  to  form 

4>r  ,  the  contribution  to  the  radial  flux  at  position  i  of  frequency  group  j: 

U 


and  thus  the  total  radiant  flux  at  position  i  s 
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Equation  (2.  9)  gives  the  expression  for  the  frequency-dependent  intensity 
to  be  used  in  Eq.  (4.  1).  The  frequency  integration  of  Eq.  (2.  9)  is  discussed 
in  Section  VII  of  Volume  V,  but  the  current  SPUTTER  code  does  not  include 
the  transmission  functions.  The  first  two  terms  of  Eq.  (2.  9)  which  form 
the  diffusion  limit  can  be  integrated,  as  in  Section  3.  3,  to  give 


i..  =  &.. 

ij  U 


in  which  the  first  term 


J£ 

<r 


R. 

J 


3B. . 
iJ 

3r 


(diffusion  limit)  , 
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is  the  frequency-group  Planck  function  and  the  second  term  contains  the 
frequency- group  Rosseland  mean  absorption  coefficient  <rg.  8  p  tj  .  In 
this  form,  Eq.  (4.2)  correctly  gives  the  frequency-group  intensity  for 
the  optically  thick  limiting  case.  The  remaining  Bj  and  aB/a-r^  terms  of 
Eq.  (2.  9)  are  formed  in  the  same  way.  Thus, 


(4.  3) 


In  Eq.  (4.  3),  mean  values  of  the  exponentials  have  been  extracted  from  the 
frequency  integrals  and  the  outstanding  problem  is  to  specify  their  values. 
Two  options  are  available;  they  differ  in  the  absorption  coefficient  used  to 
calculate  the  optical  depth.  The  first  is 
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and  the  second  is 
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where 
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dv 


and 


6  =  x.  -  x.  . 
i  1-1 


For  small  optical  depth,  the  correct  result  melees  use  of  the  Planck  mean 
absorption  coefficient.  From  Eq.  (2.  15)  the  frequency  integration  then  gives 
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(4.  5) 


The  above  prescriptions  for  frequency- group  means  are  far  from 
satisfying  and  call  for  further  work.  Considerable  economies  can  be  made 
through  reductions  in  the  number  of  frequency  groups  if  a  more  accurate 
means  of  averaging  wi  hin  groups  can  be  found.  Presently  used  choices  of 
frequency  groups  appear  to  give  a  reasonably  accurate  result,  however,  as 
indicated  by  comparisons  between  calculations  with  the  nominal  number  of 
frequency  groups  and  calculations  with  a  very  large  number  of  frequency 
groups.  (It  is  expected  that  a  unique  correct  result  will  be  obtained  as  the 
number  of  frequency  groups  is  increased,  irrespective  of  the  choice  of  the 
weighting  function  in  the  frequency-group-average  absorption  coefficient. ) 
Consequently,  a  very  few  frequency  groups  should  be  adequate  if  a  suitable 
averaging  procedure  were  developed. 

Even  with  a  crude  averaging  scheme,  considerable  improvement  in 
accuracy  results  from  choice  of  frequency-group  boundaries  so  as  to  reduce 
the  variation  of  the  absorption  coefficient  within  the  group. 

Work  on  the  absorption  coefficient  for  air  indicates  that  approximately 
20  groups,  carefully  selected  as  to  their  locations,  afford  quite  adequate 
resolution.  Enough  information  is  known  about  air  to  make  this  selection 
appear  quite  reasonable.  Air  absorption  coefficient  tapes  (DIANE)’11  have 
been  prepared  for  18,  20,  and  90  groups.  The  90~group  tape  is  used  to 
check  on  the  frequency  integrations  at  selected  times.  The  proper  averages 
to  use  are  difficult  to  decide  on  at  this  time.  There  are  provisions  for 
reading  into  storage  from  the  DIANE  tapes  both  the  Rosseland  and  Planck 
averages,  which  are  used  at  present  in  the  thick  oi  thin  limits,  respectively. 

$ 

See  Section  VI  of  Volume  V. 
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SECTION  V 

SUBROUTINE  ORGANIZATION  AND  ECONOMICS 


The  present  spherical  transport  subroutines  were  written  with  the  idea 
of  removing  unnecessary  calculations  from  inside  the  frequency  loop  and 
characteristic  ray  (y-line)  integrations  while  improving  the  accuracy  in  the 
angular  integrations  by  using  an  interpolation  between  y-lines.  These 
improvements  required  an  increase  in  storage  for  the  subroutines  to  attain 
a  decrease  in  calcuiational  time  and  an  increase  in  accuracy.  It  is  now 
practical  to  use  fewer  y-lines  and  thus  a  factor  of  approximately  4  in  savings 
on  calcuiational  times  over  the  old  routines  may  be  achieved.  The  reorganized 
subroutines  will  be  discussed  in  two  sections,  corresponding  to  the  two 
major  subroutines:  (1)  the  radiation  subroutine  (SRADTN)  in  which  most  of 
the  preliminary  setup  and  the  diffusion  calculation  is  completed  and  (2)  the 
transport  subroutine  (STRANS)  in  which  the  intensity  calculation  and  angular 
integrations  are  performed.  The  subroutines  which  execute  the  opacity 
interpolations  (KAPPA),  Planck  function  (PJLNKUT^2)),  and  fast  exponential 
(FREXP)  will  be  discussed  in  Section  5.  3.  The  input  numbers  and  the  output 
edits  will  be  presented  in  Sections  5.  4  and  5.  5. 

5.  1.  THE  SRADTN  SUBROUTINE 

In  SRADTN,  the  main  y-line  array  is  set  up,  high-frequency  groups 
are  merged,  a  source  region  is  established,  boundary  sources  and  derivatives 
are  calculated,  regions  for  transport  and  diffusion  are  formed,  diffusion 
fluxes  are  calculated,  frequency  integration  is  performed,  and  the  radiation 
time-step  control  is  evaluated.  Each  of  these  activities  in  SRADTN  will  be 
discussed  in  subsequent  paragraphs. 

5.1.1.  Set  Initial  Y-line  Array 

The  number  of  y-lines  is  limited  by  the  storage  allocated  for  the 
x-block.  A  test  is  made,  and  if  this  storage  is  to  be  exceeded  by  the  placing 
of  a  y-line  at  each  radius,  the  y-array  is  reconstructed  by  using  every  other 
radius  for  the  y-line  placement.  However,  an  additional  y-line  is  added  at 
each  zone  interface  where  the  temperature  gradient  is  large  (refer  to 
Section  2.2.  1 ), 
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5.1.2.  Merge  Frequency  Groups 

Frequency  groups  that  are  too  far  out  on  the  Planck  tail  for  a 
"maximum"  temperature  in  the  mesh  are  merged.  The  criterion  used  is 
as  follows:  If  the  lower  frequency  boundary  hvj  of  the  group  in  question 
(hVj,  hv^)  is  greater  than  ten  times  the  maximum  temperature  (THMAX) 
in  the  mesh,  this  group  will  be  merged  with  the  next  lower  group.  Merging 
will  continue  until  over  half  the  groups  have  been  merged;  at  this  point, 
either  the  calculation  is  terminated  or  a  second  DIANE  tape  is  called.  On 
merging,  Rosseland  and  Planck  averages  are  formed  by  using  the  following 
equations  for  dB/d0^  and  the  appropriate  sums: 
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The  Planck  weighting  functions  {bj)  are  obtained  from  PLNKUT,  as  described 
later.  On  completing  the  merging,  the  merged  opacities  are  formed: 


^dBv/d0^/^dBv/(d04 


(CAPAR)  , 


(CAPAC)  • 


(5.2) 


5.  1.  3.  Set  Up  Sources  and  Derivatives 

The  frequency-dependent  sources  must  be  established  at  the  interfaces 
from  the  zonal  quantities  bj0^+i  (X6(i))  and  T^+i  (H3(i)).  The  difference 
equations  used  were  given  in  Section  2.  1.  Before  the  calculation  of  the 
Planck  function  (bj)  is  made,  i,  e.  ,  before  calling  PLNKUT,  a  test  is  made 
to  see  if  U|  (i.  e.  ,  Jthe  reduced  frequency  hv  j/0)  >  19,*  if  so,  b,  =  0  (i.  e.  ,  the 
source  X6(i)  =  0.  0).  If  uj  <  19.  0  and  U£  <  0.  01,  then  bj  =  0  also,  assuming 
that  for  0^  <  10^,  the  small  bj  (bj  ~  10-5)  will  produce  a  negligibly  small 
source  contribution.  An  index  (ICX)  is  set  equal  to  the  last  zone  that  contains 
a  source,  This  source  index  is  used  to  limit  the  transport  calculation  to 
the  region  containing  sources.  While  setting  up  the  sources  and  derivatives, 
tests  are  made  on  their  discontinuous  nature  to  use  either  a  linear  or  constant 
form  in  the  intensity  integrations.  The  initial  check  is  on  the  minimum  optical 
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depth  of  adjacent  zones  to  ensure  that  both  are  transparent  (less  than  0.  3). 

If  this  condition  holds  and  if  both  the  sources  and  optical  depths  are  changing 
rapidly  in  r  (change  greater  than  a  factor  of  two),  the  derivative  at  that 
interface  (TG(i))  is  set  equal  to  zero.  The  zero  source  derivative  is  used 
in  STRANS,  as  a  test,  to  set  up  the  constant  source  terms.  For  the  intensity 
integration,  special  boundary  sources  and  derivatives  are  also  established 
at  the  edge  of  the  source  region  (I  =  ICX)  and  at  the  outside  of  the  mesh 
(I  =  IM)  (see  Section  2.  1.4), 

5.1.4.  Determine  Diffusion  Region 

The  principal  criterion  for  defining  a  diffusion  region  is  that  the  first 
derivative  of  the  source  function  (TG)  be  small  compared  to  the  source  (Y2) 

(see  Section  3.2).  When  the  zone  is  found  to  be  diffusion,  the  boundary  is 
tagged  by  setting  (X3  =  -1.  ).  Before  incorporating  this  interface  into  a 
diffusion  region,  the  possible  influence  from  sources  on  either  side  is 
considered  and  a  further  test  is  made.  From  the  last  diffusion  boundary, 
a  test  is  made  for  an  optical  depth  in  succeeding  zones  to  the  left.  If  more 
than  HVB  optical  depths  appear  in  the  next  zone,  then  this  zone  is  calculated 
by  transport  and  removed  from  the  diffusion  region  (set  X4(i)  =  -1,  0).  HVB 
is  an  input  number,  which  is  usually  around  5.  When  r  =  0  is  reached  after 
testing  each  zone,  zones  out  to  the  right  of  the  present  transport  region  are 
tested  in  the  same  manner.  The  above  test  buffers  the  transport  region 
with  an  (HVB)  mean-free -path-thick  diffusion  boundary.  If  the  zone  boundary 
stays  diffusion,  i.  e. ,  X3(i)  =  -1,  0  and  X4(i)  =  0.  0,  a  diffusion  flux  is  calcu¬ 
lated  from  the  source  gradients,  as  described  in  Section  3.  1.  The  regions 
where  X3(i)  =  0.  or  X3(i)  =  -1.  and  X4(i)  =  -1.  have  been  established  as 
transport  regions  because  they  did  not  meet  the  diffusion  criteria  or  they 
reverted  to  transport  regions  by  the  optical-depth  test  described  above. 

This  transport  region  is  then  identified  by  setting  the  left  boundary  to  IAX 
and  the  right  boundary  to  IBX.  More  than  one  trans  region  may  be  set  up 
in  SRADTN,  and  if  so,  an  STRANS  calculation  will  be  made  for  each  region. 

No  one-zone  diffusion  region  is  allowed  and  the  region  outside  the  sources 
(I  >  ICX)  is  always  considered  a  transport  region. 

5.  1.  5.  Time- step  Control  and  Monofrequency  Calculation 

These  two  aspects  of  the  new  code  are  related  since  the  "grey"  absorption 
coefficients  from  the  DIANE  tape  are  used  to  estimate  a  radiation  time  step 
as  well  as  to  form  the  monofrequency  time -dependent  calculation.  In  the 
multifrequency  calculation,  after  all  groups  have  been  processed,  an  addi¬ 
tional  call  for  KAPPA  is  made  to  read  in  the  grey  absorption  coefficients. 

These  averages  were  obtained  by  integrating  the  frequency-dependent 
absorption  coefficients  for  both  Planck '(/Cp)  and  Rosseland  ( /c ^)  in  the  DTANE 
code.  The  actual  time  step  for  radiation  transfer  is  then  obtained  from  the 
formula 
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(5.3) 


At.  =  (0.  5  +  1.  5  H3(i)V(acico03)  x  CV(i)  , 

K  K. 

where  CV(i)  is  the  specific  heat  and  ac  =  4.  12  x  10*^.  The  m£  point  in 
question  is  also  checked  to  ensure  that  it  will  not  gain  or  lose  more  than 
half  its  original  energy: 

AtR  =  0.  5  x  CV(i)  x  Q(i)  x  G(i)/  |  ER(i)  \  ,  (5.  4) 

where  ER(i)  is  the  divergence  of  the  flux  and  G(i)  is  the  mass  in  the  zone. 

The  minimum  of  these  values  is  compared  to  the  hydro  time  step  (Courant) 
and  if  smaller, 

NRAD  =  F1X(DTH2/DTRMIN)  and  DTR  =  DTH2/NRAD  (5.  5) 
is  set  to  cycle  NRAD  times  through  the  radiation  routine. 

The  monofrequency  calculation  also  uses  the  grey  absorption  coefficients 
from  the  DIANE  tape.  If  KMAX  =  0.  0  and  SI  5  =  1.0,  the  frequency -averaged 
opacities  are  bypar  sed  on  the  tape  and  only  the  grey  absorption  coefficients 
are  read  into  storage.  For  succeeding  cycles,  SI 5  is  set  equal  to  zero  and 
the  interpolations  for  /c^  and  Kp  are  performed  in  KAPPA  using  the  stored 
opacities  originally  read  into  KAPPA's  common  storage.  When  the  problem 
is  restarted  it  is  therefore  necessary  to  reload  S15  equal  to  one.  If  the 
DIANE  tape  is  not  designated  (the  tape  unit  assigned  must  be  stored  in 
AMASNO(J+17),  where  J  is  the  material  number),  then  the  KAP  routine  is 
called  (KAP8  for  air)  and  used  for  the  monofrequencv  calculation. 


5.  2.  THE  STRANS  SUBROUTINE 


The  subroutine  STRANS  is  called  by  SRADTN  to  carry  out  the  intensity 
integration  between  IAX  and  IBX,  saving  various  quantities  on  the  inward 
pass  that  will  be  used  on  the  outward  pass  as  well  as  the  angular  integration 
of  the  flux  between  y-lines  (if_{  tyi  d/i).  The  logic  is  to  calculate  first  the 
intensities  along  the  central  ray  (y  =  0),  storing  these  values  and  the 
exponentials  (e"^T)  for  use  on  the  outward  pass.  The  outward  pass  is  then 
calculated  using  the  above  exponentials  and  the  new  intensities  are  s  .  ;  ed 
separately.  After  the  intensity  transport  along  a  typical  ray  (y  £  0)if,  done, 
the  flux  calculation  is  done  while  the  right  side  of  the  calculation  is  being 
completed.  The  angular  integration  is  based  on  a  linear  interpoiation  of 
the  intensities  between  y-lines.  The  logic  in  STRANS  is  described  in  detail 
in  the  following  sections. 
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5.  2. 1.  Y-line  Placement 


In  SRADTN,  a  complete  set  of  y-lines  was  established,  based  on 
storage  limitation  and  temperature  gradients.  In  STRANS,  y-lines  in  each 
frequency  group  are  selected  more  judiciously,  based  on  a  number  of  conditions. 
First,  if  a  diffusion  core  exists  (LAX  >1),  three  rays  that  will  penetrate 
the  core  are  selected.  As  each  ray  is  selected,  the  intensity  calculation 
and  angular  integration  is  completed.  In  the  transport  region,  a  test  is 
made  for  steep  gradients  between  pairs  of  zones,  as  indicated  in  Section  2.2.2. 
If  changes  greater  than  approximately  15  percent  occur,  a  y-line  is  selected; 
if  not,  every  fifth  y-line  is  used. 

5.  2.  2.  Intensity  Integration  Along  Y-lines 

As  indicated  in  Section  II,  the  central  y  =  0  integration  is  carried  out 
first  and  the  intensities  at  each  intersection  are  stored  in  working  arrays. 

The  integration  using  Eq,  (2.  9)  starts  from  the  outside  IBXP1  and  works 
inward,  storing  the  exponentials  e"^T  in  H4(i)  and  the  intensities  I.  in 
SUMX3(i)  to  zone  LAX.  For  the  outward  pass,  these  exponentials  are  used 
in  Eq.  (2.  9),  but  now  the  sign  of  the  derivatives  is  changed  to  dB/dr  -*■  -dB/dr. 
After  the  selection  of  a  y-line,  the  intensity  calculation  is  performed.  Now, 
not  only  the  exponential  is  stored,  but  also  ji(8B/8h)-+j.  in  X8(i)  on  the  inward 
pass  for  use  on  the  outward  pass.  During  the  outwarcf  pass,  the  angular 
integration  is  accomplished  and  stored  in  X2{i)  using  the  intensities  calcu¬ 
lated  on  both  y-lines,  as  described  in  Section  5.3.  The  special  condition 
where  a  y-line  intersects  a  zone  near  x  =  0  is  treated  differently  (see 
Section  2.  1). 

The  regions  where  constant  sources,  and  therefore  zero  boundary 
derivatives,  should  be  used  in  the  intensity  integrations  were  established 
in  SRADTN  by  setting  TG(i)  equal  to  zero.  In  the  integration  along  a 
particular  y-line,  a  test  is  made  on  TG(i)  ac  each  interface;  if  zero,  the 
sourqe  terms  Y2(i  -  1)  and  Y2(i)  are  set  equal  to  X6(i  +  £)  (see  Fig.  2.  2). 

As  was  discussed  in  Section  2.  1,  1,  the  accuracy  of  the  exponential 
term  and  the  effect  ox  truncation  errors  mean  that  the  general  formula  will 
not  reduce  in  the  limit  of  small  optical  depths  to  the  transparent  case.  To 
>orrect  this  situation,  a  test  is  made  on  t.  (the  one-half  optical  depth  stored 

1  y 

in  H3(i)  as  calculated)  and  if  it  is  less  than  10  ,  a  switch  is  made  to  the 

limiting  form  of  the  transport  equation  (Eq.  (2.  15))  developed  in  Section  II. 

5.  2.  3.  Angular  Integration 

This  section  of  STRANS  uses  the  intensity  information  stored  on  adjacent 
y-lines  to  give  a  more  accurate  angular  integration.  The  iinear  interpolation 
of  the  intensities  is  used  in  the  quadrature  formulation  of  the  fluxes  (Eq.  (2.32) 
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in  Section  2.  3).  While  the  intensity  integration  is  being  completed  on  the 
right-hand  side  (RHS),  use  in  made  of  stored  intensities  on  the  left-hand 
side  (LHS)  for  adjacent  y-lines  to  calculate  the  fluxes  across  radial  bounda¬ 
ries  X2(i).  The  stored  quantities  used  are:  SUMX3(i)  for  the  LHS  intensity 
on  the  present  y-line,  FM(i)  for  the  LHS  intensity  on  the  previous  y-line, 
and  SUMX4(i)  for  the  RHS  intensity  from  the  previous  y-line  pass.  F2  is 
the  current  intensity  being  calculated  along  the  y-line  in  question.  The 
special  condition  to  complete  the  flux  integrations  where  no  y-lines  exist 
uses  a  "top-slice"  approximation,  as  discussed  in  Section  5.  2.  5. 

5.  2.  4.  Boundary  Concutions 

The  usual  assumption  of  zero  intensity  at  the  surface  of  the  sphere 
is  used;  i.  e. ,  L_  ^  =  0  to  start  the  intensity  calculation.  If  the  indication 
is  that  a  diffusion  region  bounds  the  transport  region  (IBX  <  IM),  a  diffusion 
boundary  condition  for  the  intensity  is  applied  (see  Section  2.  1.  5).  At  the 
center,  the  outward  intensity  in  set  equal  to  the  inward  intensity. 

5.  2.  5.  "Top  Slices"  and  Finish  Up 

The  term  "top  slice"  is  used  to  describe  the  integration  of  the 
flux  for  zones  which  do  not  have  y-lines.  For  instance,  outside  the 
sources  there  is  no  need  for  y-lines,  and  therefore  it  is  assumed  that  the 
intensity  drops  to  zero  at  the  outside.  If  there  exist  zones  with  sources 
(the  normal  case)  between  y-lines,  then  a  linear  interpolation  from  y^  to 
Y2  is  used  to  establish  intensities  for  integration  at  those  zones.  The  inter¬ 
polation  form  used  is  given  in  Section  III.  If  a  diffusion  region  bounds  the 
transport  region,  a  diffusion  intensity  for  the  interpolation  is  calculated 
(see  Section  III). 


5.  3.  AUXILIARY  SUBROUTINES 


In  addition  to  the  two  new  basic  subroutines  SRADTN  and  STRANS, 
some  changes  have  been  made  in  the  auxiliary  subroutines  EXP,  PLNKUT, 
and  KAPPA.  These  changes  include  (1)  a  fast  exponential  (FREXP),  (2)  a 
two-argument  Planck  function,  and  (3)  the  use  of  the  average  opacities  from 
KaPPA  (0  and  p  interpolation)  for  the  monofrequency  calculation  as  well  as 
for  the  Planck  opacities. 

The  new  fast  exponential  routine  uses  table  lookup  and  interpolation 
rather  than  the  normal  expansion  methods.  The  routine  is  written  in 
machine  language  but  uses  the  library  routine  EXP(X)  for  positive  X  or 
X  >  -10.  An  over-all  gain  in  speed  of  a  few  percent  was  achieved  in  one 
comparison  SPUTTER  calculation. 
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The  PLNKUT  routine,  with  its  associated  tables  PLNKTT,  has  been 
corrected  and  made  more  efficient  by  using  a  two-argument  call  which  now 
calculates  from  either  the  analytic  form  or  from  the  tables  the  difference  in 

^  1  f  V2  hv3  1  j  ,r  ,, 

b<ullU2)  =  B|  —  -WkT  ;dl,sbj  •  (5-6) 

JVj  c  e  -  1  J 

The  accuracy  is  improved  since  now  not  only  differences  of  nearly  equal 
numbers  are  subtracted. 

The  subroutine  KAPPA,  which  calls  in  the  group-averaged  abeorption 
coefficients  from  the  DIANE  tape  and  performs  a  bilinear  log  interpolation 
in  temperature  and  density,  has  been  modified  to  obtain  the  grey  absorption 
coefficients  as  well  as  the  Planck  averages.  At  present  the  format  of  the 
DIANE  (absorption  coefficient)  *ape  includes  a  BCD.  record  for  tape  identi¬ 
fication,  the  Rosseland  and  Planck  averages  for  a  selected  set  of  tem¬ 
peratures  and  densities  from  0.  25  ev  to  50  ev  and  from  10  normal  to  10"^ 
normal,  and  the  actual  integration,  J  /cy  dv,  for  the  grey  case.  The  grey 
or  frequency-integrated  averages  are  also  used  for  an  estimate  of  the  time 
step  in  SRADTN.  KAPPA  reads  in  first  the  tape  name,  the  number  of 
frequency  groups,  and  the  size  of  the  records.  If  the  sentinel  for  multi¬ 
frequency  is  set  to  KMAX  =  1,  then  the  first  frequency  group,  hvj,  and  its 
absorption  coefficients  are  read  into  storage.  The  interpolations  in  log  0^ 
and  log  pj  are  performed  and  a  return  to  SRADTN  is  made.  If  KMAX  =  0, 
then  KAPPA  skips  over  the  frequency-dependent  absorption  coefficients 
and  reads  into  storage  the  grey  averages.  A  signal,  SI 5  =  0,  is  subsequently 
set,  and  for  further  cycles  the  interpolations  are  made  on  the  stored  quantities; 
th'e  tape  is  not  called  again. 


5.4.  INPUT  NUMBERS 

The  input  quantities  used  in  the  radiation  transport  subroutines  and 
their  functions  are  listed  below  by  card  number.  Included  is  a  set  of  values 


for  the 

input  quantities  selected  for  solving  typical  problems. 

Card 

No. 

Quantity 

Value 

Purpose 

44 

KMAX 

1. 

Signal  to  do  multifrequency 

77 

CVB 

0 

Select  all  y-lines 

81 

HVB 

5. 

Buffer  of  trans  region  in  number  of  optical  depths 

83 

HCB 

0.  1 

Diffusion  criteria  (see  Section  5.  4) 
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Card 


No. 

Quantity 

Value 

Purpose 

87 

CB 

0 

Brightness  print  on  cycles 

88 

GA 

0.  33 

Gradients  in  optical  depths  and  sources  for 
applying  TG  criteria  (see  Section  5.  3) 

90 

GL 

0.3 

Gradients  in  source  and  MFP  for  selecting 
y-lines 

121 

CMIN 

0.3 

Minimum  depth  of  adjacent  zones  for  applying 
TG  criteria  (see  Section  5.  3) 

127 

AC03T4 

0.  1 

One-half  optical  depth  for  Planck- Rosseland 
transition 

147 

S12 

10. 

Cycles  between  multifrequency  prints 

150 

S15 

1. 

Restart  on  grey  calculation 

8478 

TELM(37) 

0.005 

Fraction  of  total  energy  in  zone  for  time- step 
criteria 

8725 

AMASNO(25) 

9. 

To  select  B5,  in  this  case  for  air  tape 

8858 

SOLID(IO) 

10. 

Signal  to  use  Planck  mean  free  paths 

8466 

TELM(25) 

1. 

Constant  multiplying  radiation  time  step 

5.  5. 

EDITS 

The  editing  of  such  frequency-dependent  quantities  as  H3,  the  optical 
depth  (Rosseland),  X6,  the  source  (b.0^),  X2,  the  flux  {in  ergs/4/3  ir  sec), 

ER,  the  radiation  energy  (in  ergs/crfP),  PR,  the  radiation  pressure  (in 
dynes/ cm^),  PHO,  the  sidewise  flux  (in  ergs/ cm^- sec),  PK1,  the  forward 
flux  (ergs/cm2-sec),  and  PH2,  the  backward  flux  (in  ergs/ cm^-sec)  versus 
radius  is  accomplished  by  setting  S12  to  the  desired  number  of  cycles  between 
prints.  These  multifrequency  edits  have  been  used  to  evaluate  the  criteria 
for  the  subroutines  as  well  as  for  diagnostics  during  the  calculations.  For 
example,  one  quantity  found  to  be  moot  useful  for  fireball  diagnostics  has 
been  the  flux  out  in  various  frequency  g?  oups.  Because  of  this,  the  flux  out 
is  edited  out  in  the  standard  SPUTTER  output,  X7. 

\  list  of  sample  editing  for  a  particular  frequency  group  is  given  on 
page  45.  The  HNU  is  in  electron  volts. 

To  obtain  a  limited  set  of  intensities  for  purposes  of  debugging  the 
code,  a  "debug  print"  option  can  be  added.  The  FORTRAN  statements  and 
input  cards  necessary  for  debugging  are  listed  on  page  46.  The  quantities 
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found  useful  to  display  for  each  frequency  group  and  for  a  central  y-line 
and  a  selected  y-line  are  listed  on  pages  47  and  48. 
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or»r»  non  r»r> 


Table  5.  2 

INPUT  CARDS  FOR  DEBUG  OPTION 


* 


♦ 


LOCATION  OF  VARIABLE  IN  SPLITTER  COMMON 
NUMBER  OF  INPUT  WORDS  THIS  CARO 
♦  INPUT  WORD 


311  1  0. 

6876  1  5.9 

8877  i  6.1 


R! 1521 
SOLID! 28 J 
SOLID! 29) 


DEBUG  PRINT 

70S  DHNU*HNUP—HNU 

IF!R(152))713, 712,713 

713  IF  f HNU— SOLID! 291 )  712,714,714 

714  IF  ! HNU- SOLID ! 28 ) )  715,715,712 

715  WRITE  16,31  HNU,HNUP 

717  WRITE  f 6,5)  Y(JJ) 

718  WRITE  !6,6)  JJ, IAX, IBX, ICX 
WRITE  !6,8) 

IAXP1*MAX0< IAX»IAXP—i) 

DO  721  I*IAXP1, IBXP1 
IF! JJ.GT.l)  GO  TO  720 

PRINT  Y*0  INTEGRATION 

719  WRITE  !6,7)  I,C!I),X6!I) ,H2! I) ,H3! I) ,SUMX3( I ) , SUMX4i I ) ,X2( I ) 
GO  TO  721 

720  KXK-KK-I4IM43 

PRINT  REGULAR  Y-LINE  INTEGRATION 

WRITE  16, 7)1,  X(  KXK),  X61I), H2!  1 1  ,H3(  I  ),SUMX3!  I)  »SUMX4!  I), X2U) 

721  CONTINUE 
712  CONTINUE 
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SECTION  VI 

TIMING  AND  ACCURACY  COMPARISONS 


A  timing  and  accuracy  study  was  made  on  the  new  radiation  and  trans¬ 
port  subroutines,  SRADTN  and  STRANS,  so  that  the  new  routines  could  be 
compared  with  the  old  routines  to  obtain  a  quantitative  estimate  of  the 
improvements.  It  was  estimated  that  the  reorganization  of  the  subroutines 
would  reduce  the  calculational  time  by  50  percent. 


6.  1.  CALCULATIONS  FOR  TIMING  AND  ACCURACY 

1.  Initial  formulation  of  sources,  including  source  at  icx  +  1,  and 

characteristic  y-lines  was  made  with  the  faster  exponential 
subroutine  FREXP.  (This  is  compared  with  the  same  calculation 
using  the  library  exponential  subroutine.  )  I 

2.  In  the  subsequent  calculations,  the  FREXP  was  used  and  sources 
at  icx  +•  1  were  set  to  0. 

3.  The  index  for  setting  the  initial  y-line  distribution  (all  y-lines 
are  not  necessarily  used)  was  improved.  A  calculation  was  made 
where  the  input  number,  CVB,  that  controls  the  initial  number  of 
y-lines  was  set  to  0  (i.  e.  .  generate  all  y-lines). 

4.  The  above  calculations  were  done  with  thick/ thick  (i.  e.  ,  Ross  eland) 
opacities  used  throughout  (SOLID(IO)  =  10.  ).  A  two-cycle  run  was 
completed  to  check  that  the  calculation  recycles  correctly. 

5.  The  main  accuracy  and  timing  comparisons  were  made  with  cal¬ 
culations  having  normal  y-line  and  all-y-line  configurations  and 
with  the  thick/thick  and  thick/thin  approximations  (set  of  4).  The 
thick/ thin  approximation  used  here  differs  from  the  Planck  and 
Roaseland  averages  used  in  the  old  subroutines. 

The  comparison  in  calculation  time  was  obtained  by  using  timing  calls 
at  selected  locations  in  the  logic  of  the  code.  To  use  the  timing  calls,  it  was 
necessary  to  establish  a  fiducial  time  from  the  system  clock  and  then  print 
the  location  of  the  time  call,  the  time, .  and  the  difference  in  time  between 
calls  for  each  calL  The  subroutine  that  carries  out  these  steps  is  CLOCK. 

In  the  calculations  described  above,  the  subroutine  CLOCK  was  called 
at  the  following  locations  in  SRADTN  and  STRANS: 
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SRADTN 


13.  105  -  After  setting  logjg  of  9  and  SV 

13. 118  -  After  calculating  x-array 

13.  140  -  After  call  KAPPA  on  merge 

13.  700  -  End  of  merge  and  source 

13.  701  -  After  call  KAPPA  on  main  frequency  loop 

13.  151  -  After  calculating  general  sources 

13.  180  -  Before  call  transport  (STRANS) 

13.292  -  After  edit  (normal)  end  of  frequency  loop 
13.286  ~  After  last  frequency  start  time  step 

13.  289  -  End  of  cycle 

STRANS 

14.  708  -  Before  debug  print 

6.  2.  COMPARISON  FOR  TIMING 


The  general  conclusions  from  the  timing  study  are  summarized  as 
follows: 

1.  The  use  of  the  faster  exponential  subroutine  FREXP  in  comparing 
a  normal  y-line  problem  showed  an  improvement  of  approximately 

2.  5  percent. 

2.  The  all-y-line  calculations  did  not  increase  the  calculational  time 
in  proportion  to  the  increase  in  y-lines,  i.  e.  ,  a  factor  of  5  increase 
in  y-lines  increased  the  running  time  by  only  80  percent. 

3.  The  normal  setup  times  and  typical  frequency-group-dependent 
calculational  times  for  a  number  of  runs  (in  units  of  1/60  sec) 
are: 


Set  up  logjQ  Q  and  SV  ~4 

Calculate  x-array  (36  y-lines)  ~32 

Call  KAPPA  first  group  ~34 

Call  KAPPA  second  group  ~23 

Merge  above  and  calculate  source  ~10 

Total  setup  time,  including  merging 
two  groups  ~103 
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4.  A  breakdown  of  the  time*  required  for  a  normal  frequency  loop 
(i.  c. .  nc-rmal  y-lines)  is  as  follows: 

Bring  in  new  absorption  coefficients 


(call  KAPPA)  ~18 

Source  calculation  (PLNKUTT)  ~5  -  8 

Calculations  before  "call  trans.  " 
mostly  diffusion  fluxes  ~2 

y-line  calculation  of  intensities  ~7 

Normal  print  {X2,  H3,  etc. )  ~35 

Total  time-  per  frequency  group 
with  normal  EDIT  ~68 

Brightness  print  ~11 


The  radiation  cycle  is  completed  by  calculating  a  minimum  time  step 
from  the  grey  absorption  coefficients  obtained  at  this  time  from  a  fit  in 
KAP8. 

The  total  times  for  a  radiation  cycle  with  ~60  zones  (36  are  active) 
and  21  frequency  groups  require  for 

3 

Normal  editing  -*■'1.  49  x  10  (1/60  sec) 

Brightness  editing  ~i.  2  x  10^  (1/60) 

No  special  editing  ~1.  00  x  10^  (1/60) 

3 

The  ali-y-line  calculation  with  normal  editing  takes  ~2.  40  x  10  (1/60  sec). 

In  comparison,  with  the  old  subroutines  the  calculation  times  (in 
units  of  1/60  sec)  for  a  particular  frequency  group  require: 


All  y-line 

~200 

KAPPA 

~17 

Source 

~1 3- 30 

Normal  prints 

~1 7 

Brightness  print 

~13 

The  setup  time  for  all  groups,  including  the  monofrequency  time  step,  is 
~30  (1/60  sec).  The  total  times  for  a  radiation  cycle  with  60  zones,  21 
frequency  groups,  and  all  y-lines  require  for 

3 

All  prints  ~5.  40  x  10  (1/60  sec) 

Normal  editing  ~5.  19  x  10^  (1/60  sec) 
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6.  X  COMPARISON  FOR  ACCURACY 


The  new  SRADTN  and  STRANS  subroutine  calculations  are  compared 
to  an  all-y-line,  all-transport  result  from  HA5.  In  general,  the  results 
all  agree  within  5  to  S  percent  either  from  a  flux  versus  radius  comparison 
(Fig.  6.  1)  or  from  a  central  brightness  comparison  (Fig.  6.2).  The  dif¬ 
ferences  are  qualitatively  understood  and  are  within  acceptable  limits.  The 
close  agreement  (within  1  percent)  of  the  new  all-trans,  all-y-line  calcula¬ 
tions  with  the  standard  routines  indicates  that  the  normal  y-line  calculation 
for  Rosseland  averages  (~5  to  8  percent)  gives  very  gooJ  results  with 
fewer  y-lines.  This  result  can  probably  be  improved. 


6. 4.  CONCLUSIONS 


The  over-all  computational  time  has  been  reduced  by  a  factor  of  two 
from  the  old  standard  transport  calculation  without  apparent  loss  in  accuracy. 
The  new  formulation  will  allow  for  fewer  y-lines  and  subsequent  further 
decrease  in  running  time  from  the  old  calculation  with  estimated  differences 
of  only  5  to  8  percent. 
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AIBFTC  SRADTN  Li ST#06CK#R£>r 
SUBROUTINE  RAOTN 
C  DATE  COMPILED  3  JUNE  1965 

C  SPECIAL  EDIT  ON  S12 

C  REVISED  SPHERICAL  TRANSPORT  (NEW  65) 

C  MUST  USE  WITH  LATEST  KAPPA  FOR  AVERAGE  OPACITIES 
C  H2(I)  NOW  CONTAINS  PLANCK  MFP  ON  SOL IO( 10 )  =0*0 

C *********************************************************************** 
C*  ** 


S  P 

U  T  f  E 

R  C  0  M  M  0 

N 

♦* 

COMMON 

LM0A(37>, 

NR  , 

NSMLR  ,  IA 

,  IB 

,  ICA  , 

♦ 

ICB  , 

1 

KMAX 

,  8LANK1, 

8LANK2. 

8LANK3,  IAP1 

,  IBP1 

,  ICAP1  , 

ICBP1  , 

2 

II 

»  IG 

NRAD  , 

BLANK4,  I AMI 

,  IBMl 

,  ICAM1  , 

ICBMl  , 

3 

IIP! 

,  IGMI  , 

I  ALPHA, 

8LANK5,  TH 

,  TMAX 

,  BLANK6 , 

DELPRT, 

4 

FREQ 

,  CNTMAX, 

AR  , 

ASMLR  ,  PUSHA 

,  PUSHB  ,  BOILA  , 

BOILB  , 

5 

CVA 

,  CVB 

t 

SLUG  , 

ALPHA  ,  HVA 

,  HVB 

,  HCA  , 

HC3  , 

6 

EMINA  ,  EMIMB  , 

CA  , 

CB  ,  GA 

,  GB 

,  GL  , 

GR  , 

7 

RHOL 

,  RKCR  , 

EPIO  , 

EPSI  ,  RIA 

,  RIB 

,  RDI A  , 

RDIB  , 

8 

RPIA 

,  RPIB  , 

RPDIA  , 

RPDZB  ,  T PR I NT 

,  TA 

,  TB  , 

TC 

COMMON 

TO 

9 

TE  , 

0TH2  ,  DTH2P 

,  DTH1 

,  DTRMIN, 

OTMAX  , 

1 

DTMAXi,  0TMAX2, 

0TMAX3, 

DTR  ,  SWITCH 

,  CO 

,  CMIN  , 

DELTA  , 

2 

GAMA 

,  WCRIT  , 

SIGMAQ, 

AC  ,  AC03T4 

,  CNVRT  ,  SUMRA  , 

SUMRB  , 

3 

ROI A 

,  ROIAMl, 

ROIB  , 

ROIBPl,  GMS 

,  SI 

,  S2  , 

S3  , 

4 

S4 

,  S5 

» 

S6  , 

S7  ,  S8 

,  S9 

*  S10  , 

Sll  , 

5 

512 

,  S13 

, 

S14  , 

S15  ,  S16 

,  S17 

,  S18  , 

S19  , 

6 

S20 

,  EO 

FO  , 

TAU  ,  ZERO 

,  R 

(152),  DELTAR(152), 

7 

AS  Q 

(152), 

RO 

(152), 

VD  (1521, 

ROD 

(152),  SMLR 

(152), 

8 

DELR 

(  37), 

P 

(152), 

PI  (152), 

P8 

(152),  PB1 

(152) 

COMMON 

P2 

(152), 

SV  (152), 

RHO 

(152),  THETA  (152), 

1 

W 

(152), 

E 

(152)  , 

El  (152), 

EK 

(152),  A 

(152), 

2 

V 

(152), 

G 

(152)  , 

0  (152), 

C 

(152),  X2 

(152), 

3 

X} 

(152), 

X4 

(152)  , 

X5  (152), 

X6 

(152),  X7 

(152), 

4 

SMLA 

(152), 

SMLB 

(152)  , 

SMLC  (152), 

SMLD 

(152),  SMLE 

(152), 

5 

EC 

(152), 

ER 

(152) , 

SMLQ  (152), 

SMLH 

(152),  BIGA 

(152), 

6 

BIGB 

(152), 

CV 

( 1 52 )  , 

BC  (152), 

BR 

(152)  ,  CHIC 

(152), 

7 

CHIR 

(152), 

CAP AC  (152), 

CAPAR  (152)  , 

CRTC 

(152),  CRTR 

(152), 

8 

CRTPC 

(152), 

GOf-R 

(152), 

FEW  (152), 

CAR 

(152),  OKLM 

(  37) 

COMMON 

TELM 

(  37), 

EKLM  (  37), 

ELM 

(  37),  FCLM 

(  37), 

1 

FRLM 

(  37), 

WLM 

(  37), 

QLM  (  37), 

AMASNOI  37),  CHRNG 

(  37), 

2 

ZPl 

(  37), 

ZP2 

(  37), 

SOLID  (  37), 

ECHCK 

(  37),  RK 

(104), 

3 

RL 

(  37), 

RHOK 

(  104)  , 

RDK  (104), 

THETAK( 104) ,  TEMP 

(  16), 

4 

HEAD 

(  12), 

MAXL 

, 

MAXLM 

c*  ** 

C  «^***«***  **  t^**********  *****t>y  ********************************  ********* 

C  * 

COMMON  /LINDLY/  HNU,SGNL  ,  IHNU  ,NHNU,  HNUP,  NT , Jj^,  IN,  DHNU, THICK, NY 
C 

COMMON  /CNTRL/  SCYCLE,  JMULT 


gWfJCS J3  PA*  WAS  BLAJK  'BKHPCBt  WAS  HOT  m«D 


C 


COMMON  /DAVIS/  FEX  <152),  ICX,  IR,  PRC152), 
X  TG<  152)  *  XC2400),  GXU52) 


C 


DIMENSION 

CSQO 

(1) 

*  £RR< 1)  , 

FM 

<1), 

H 

m* 

H2 

III 

1 

H3 

m. 

H4 

<1) 

,  PRR 

<11* 

01 

ui. 

Q2 

<u» 

03 

<11 

2 

Q37 

<ii* 

Q38 

<11 

,  SUMX2 

<11* 

SUMX3 

in* 

SUMX4 

<ii* 

X8 

<1) 

3 

V 

Y2 

i  1) 

,  YSQD 

m* 

PHO 

<d* 

PHI 

<i) 

c 


EQUIVALENCE 

<  AC03T4 

•TRDBGi. 

<8  C 

»SUMX4) , 

(8IGA 

*V 

) 

1 

(6IG8 

*  H 

)  * 

IBR 

,H3 

}* 

(CAR 

*Q37 

)• 

(CHIC 

•SUMX3) 

2 

(CHIR 

*Q  38 

If 

(CR.TR 

, SUMX2 I « 

<  EC 

*Q1 

>. 

(GOFR 

f  03 

I 

3 

<RDO 

,  ERR 

If 

<  S12  * 

EDITMFI, 

<  SMLA 

*H4 

I* 

(SML8 

*FM 

) 

4 

tSMLC 

«CSQO 

)* 

ISMLD 

*Y2 

1* 

<  SMLE 

,H2 

If 

(SMLH 

fX8 

) 

5 

ISMLR 

t  PRR 

1* 

<V 

*Q2 

1* 

<X7 

,YSQD 

)  f 

(CRTC 

,PH1 

) 

6  <CRTPC,PHO  I 


C  * 

£***************♦**♦**♦♦♦*****♦*******♦*♦**♦***********# *******♦*****!*** 


C 

OX  CONTAINS  X  FROM  THE  PREVIOUS 

Y  LINE 

c 

* 

c 

CSQD 

SAME 

AS 

SMLC 

* 

c 

EDI IMF 

SAME 

AS 

s:2 

* 

c 

ERR 

SAME 

AS 

ROD 

* 

c 

FM 

SAME 

AS 

SML8 

* 

c 

H 

SAME 

AS 

SIGB 

* 

c 

H2 

SAME 

AS 

SMLE 

♦ 

c 

H3 

SAME 

AS 

BR 

* 

c 

H4 

SAME 

AS 

SMLA 

* 

c 

PHO 

SAME 

AS 

CRT  PC 

♦ 

c 

PHI 

SAME 

AS 

CRTC 

* 

c 

PRR 

SAME 

AS 

SMLR 

* 

c 

Q1 

SAME 

AS 

EC 

* 

c 

Q2 

SAME 

AS 

V 

c 

Q3 

SAME 

AS 

GOFR 

c 

Q37 

SAME 

AS 

CAR 

* 

c 

Q38 

SAME 

AS 

CHIR 

* 

c 

SUMX2 

SAME 

AS 

CRTR 

♦ 

c 

SUMX3 

SAME 

AS 

CHIC 

* 

c 

SUMX4 

SAME 

AS 

BC 

* 

c 

TRDBG 

SAME 

AS 

AC03T4 

* 

c 

X8 

SAME 

AS 

SMLH 

* 

c 

Y 

SAME 

AS 

BIGA 

* 

c 

Y2 

SAME 

AS 

SMLO 

* 

c 

YSQO 

SAME 

AS 

X7 

* 

c  * 

C*********************************************************************** 


^*#t***************#*4M<**************************#*#**************4t*l**** 

C  * 

c  PHO  SIDE  FLUX  { ERGS/CM**2~SEC)  * 

C  PHI  FORWARD  FLUX  < ERGS/CM**2-S£C)  * 


SIDE  FLUX 
FORWARD  FLUX 


<ERGS/CM**2~SEC) 

<ERGS/CM**2-S£C) 


f 


C  PH2  BACKWARD  FLUX  { ERGS/CH**2-SECJ  * 

C  PRR  RADIATION  PRESSURE  l£RGS/CN**3J  * 

C  RHG  RADIATION  ENERGY  l ERGS I CM»*3 I  * 

C  » 

c*********************************************************************** 
NT IMES=BOILB 
IMP1*IB 
IM=I8H1 
IN=IA 
INM1*IN-1 
CALL  DYCHKIKOOFX) 

GO  TO ( 17,171  , KOOFX 
17  IR*IN 

THTAMX*.025 

104  00  105  I  *  I N  ,  I M 
X3CI)«0. 

X4!I)«0. 

X5(l)*0. 

X6UJ«0. 

X7(I )*0. 

C*******************************************************************1**** 

c  * 

C  SET  OP  FOR  KAPPA  INTERPOLATION  * 

C  * 

c******^**********#**********************************^***************** 
Ql(I)=TH£rA(I)**4 
Q37U)*AL0G(THETA(I1I 
Q38I I i *AL0G4 SVC  II I 
CSQD(I)=CIII**2 
C 

c  find  ir  largest  index  of  zone  with  theta 

C  GREATER  THAN  0.05 

c 

IF1THETA< I l-THTAMXI 100,100,101 
101  THTAMX*TH£TA(I ) 

100  IFITHETACIi-0.051 105, 105,106 

106  IR*I 

105  CONTINUE 
CSQD(IB)sC<IBI**2 

C 

C  PROBLEM  STOPS  IF  MAX  TEMP  LESS  THAN  0.05 

C 

IFI THTAMX-0.05) 108,108,107 
108  Sl=13« 108 
CALL  UNCLE 

107  K*IFIXC  < FLOAT (  IB)+2«0)/  ( 2.4E3/ FLOAT!  IR)<HFLGAT!  IRI  +  1.0I/2.I ) 

C 

C  DETERMINE  WHETHER  TO  SKIP  ZONES  WITH  Y-LINES 


111  K1=K 


mi*o.Q 
YSQDt 11=0.0 
IRN1*IR— 1 


JK*1 

€♦»♦♦♦♦♦♦♦♦♦♦**♦♦****»*♦♦»♦******♦*♦♦♦♦*♦**»»♦*♦*****♦*♦♦*♦*♦*♦♦***♦♦*** 
C  * 


C  SET  UP  Y  LINES  * 

C  ♦ 

€»♦♦♦♦♦***»**♦**♦♦♦*»****♦♦******♦♦♦*♦♦*♦*»♦*♦**♦♦♦♦*****♦*♦*♦*****♦**** 

c 

C  DRAW  Y  LINES - ONE  FOR  EACH  LARGE  TEMPERATURE  CHANGE 

C 

C  NO  Y  GREATER  THAN  CUR) 

C 

00  112  I=IN,IRK1 

IF  UBS(Q37U«-1I-Q37II))~CVB1  113*113,114 

114  JK*JK+1 

Y(  JK)=C( 1*1) 

YSQ0CJK)=CSQD(I+1) 

KX—K 

GO  TO  112 

113  IFIK1U14,  114,115 

115  K1=K1-1 
112  CONTINUE 

NY=JK 
GO  TO  116 
117  K=K+1 

GO  TO  111 

(;»♦♦»♦♦♦♦♦♦♦»♦»*♦♦♦*♦♦»*******♦♦♦♦***♦*»**♦*♦♦**♦*****♦***♦********♦***♦ 
C  ♦ 


C 

c 


FINO  COMPLETE  SET  OF  X  VALUES 
(INTERSECTIONS  OF  RADII  WITH  Y  LINES) 


* 


C  * 
C*******************************************************1********** ****** 
C 


C 

c 


FORMAT=-YSQD, -NUMBER  OF  INTERSECT IGNStX'S. 


116  K2=l 

DO  118  J=2VJK 
I=IMP1 

X(K2I=-YSQ0( Jl 

K2*K2+2 

KK=  1 

120  TS1=CSQ0( I )-YSQD<  J ) 
IF{TSl)122, 122,119 
119  X(K2I=SQRT (TS1) 

K2=K2+1 

IF (K2-2400) 1191, 1191, 117 
1191  1*1-1 

KK=KK+1 
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GO  TO  120 
122  KKK*K2-KK 

XIKKK»*-4KK-1) 

118  CONTINUE 
C 

C  FINISH  X-8LOCK  WITH  NEG  NUMBER 

C 

X(K2J=-1.0 

£**»** *******»*»*»***»»**»***»»**********»»***»»****»»»**»»*»***»**»»»»» 

c  * 

C  SETUP  FREQUENCY.  LOOP  * 

C  * 

C**************************************** ******************************* 

HNUP»4000« 

C 

C  SET  UP  MAX  FREQ  BOUNDARY 

C 

HNUP4=2.5fcE14 

IHNU=»1 

00  121  l*IN, IMP1 
ERR( I lu0*0 
PRRI  n*o.o  > 

121  SUMX2I I )=0»0 

q** ********************************************************************* 

c  ♦ 

C  TEST  FOR  MONOFREQUENCY  * 

C  * 

C ************* ********************************************************** 

IF  (KMAX  «EQ«  0  )  GO  TO  150 
140  CALL  KAPPA! INt IM) 

HNU4=HNU**4 

C*********************************************************************** 

c  * 

C  MERGE  GROUPS  HAVING  HNU  LESS  * 

C  THAN  TEN  TIMES  LARGEST  THETA  * 

C  * 

C***************** ********************************************* ********* 
IF( THTAMX— 0«1*HNUI  142,141,141 
C 

C  REJECT  TAPE  IF  MORE  THAN  HALF  OF  GROUPS  MERGE 

c 

142  IF! IHNU+IHNU-NHNUI 143,144,144 
144  Sl*13.0144 
C 

C  USE  STATEMENT  AS  IN  PLANE  CODE  NOW 

C 

C  WILL  HOPE  TO  CALL  NEXT  TAPE  IN  FUTURE 

C 


CALL  UNCLE 
143  DO  145  I=IN,IM 


8ETA=HNU/TH£TA( I) 

8£TAP=HNUP/TH£TA{I) 

DFB=PLNKUTIBETA,BETAP) 

IFIOFB  .EQ.  O.J  GO  TO  145 
TEMPI1I=DFB*Q1C I ) 

EM81=EXP(*"8ETAi 

EMB2=EXPI-BETAPI 

TEMP ( 2 ) =DFB+0 .0  384974/Q 1 < I J*C HNU4/ ( l . O-EMB i I 
1*£MB1-HNUP4/(1.0-EMB2)*EMB2) 

C 

C  FORM  NUMERATORS  ANO  DENOMINATORS  OF  MERGED  KAPPAS 

C 

X6I i )=X6( I J+TEMPI1) 

X4(  I ) =X4( I I+TEMPI 21 
X5I I l=X5< I )+CAPAC( I)*TEMP< II 
X3I I )=X3( I l+TEMPI 2I/CAPARI I I 
145  CONTINUE 
HNUP=HNU 
IHNU=IHNU+1 
HNUP4=HNU4 

IF! THTAMX-O* 1*HNU)  140,703,703 

C*********************************************************************** 

c  * 

C  MONOFREQUENCY  CALCULATION  * 

C  * 

C**#^*^****************************'************************** ********* 

150  NHNU=1 
C 

C  CALLS  KAPPA  ONE  TIME  ON  RESTART  <S15=1.0>  FOR 

C  AVERAGE  OPACITIES  FROM  DIANNE  TAPE  AND  STORED 

C  IN  KAPPA  FOR  FUTURE  INTERPOLATIONS 

C 

CALL  KAPPAi IN, IM J 
DO  704  I=IN,IM 
X5< 1 1  =  1 o 
704  X6( I >-Ql( I ) 

ICX=IR 
GO  TO  139 

141  IF  C IHNU-ll  153,702,143 
C 

C  FORM  MERGED  KAPPAS 

C 

703  DO  700  I=IN, IM 

IF  ( X6  < I ) )  210,700,200 
210  Sl=13.210 
CALL  UNCLE 

200  CAPARI I ) =X4< I I/X31 1 1 
215  CAPACII)=X5U»/X6m 
700  CONTINUE 
HNUP=4000  * 
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HNUP4*2»56£14 
IHNU*IHNlr-l 
GO  TO  702 

C*********************************************************************** 

C  ♦ 

C  TYPICAL  GROUP  CALCULATIONS  OP  SOURCES  * 

C  * 

(;«♦$*#*♦*♦*♦*♦**♦♦**♦*****♦*♦******♦♦******♦*«♦*«*♦♦♦♦***♦*♦*♦♦***♦****** 

701  CALL  KAPPAUNfIMJ 
HNU4»HNU**4 

702  ICX»IN 

00  151  ISIN» 1R 
B£TA*HNU/TH£TA( II 
EMB1*EXP<-BETA) 

C 

C  AVOID  CALCULATION  OF  DFB  LESS  THAN  IE-5 

C 

IF! 8ETA-19. 1146,147,147 

146  B£TAP*HNUP/THETA( I ) 

EHB2*EXP(-8ETAP) 

IP I BETAP-.01 J 147, 147,148 

147  IF  I  ICX~*IR)  149,157,157 
157  ICX=J-l 

149  X6( I )=0.0 


GO  TO  151 


X5(IJ=0.0 


FORM  SOURCES  X6  ANO  X5 


148  OFBsPLNKUTI BETA, BET AP) 

X6< I )*0FB*Q1( I ) 

132  ICX*IR 
151  CONTINUE 
139  Q31=0.0 

£***♦****»*********♦*******************♦***********************♦******** 

c  * 

C  FORM  ROSS  ANO  PLANCK  OPTICAL  DEPTHS  * 

C  * 

£*****«♦♦«*♦♦**♦♦**♦******♦**♦*♦♦***♦***♦♦♦*♦****♦«****♦♦**♦♦*****♦**♦♦* 
00  152  I*IN,IM 
IF(CAPAR( I ) ) 153,153,154 
154  IF4CAPAC( 11)153, 153, 155 
153  Sl*13« 153 
CALL  UNCLE 


CHOOSE  ALL  ROSS  IF  SOLID  10  IS  POS. 


155  IF (SOLID! 10) «EQ«0« )G0  TO  167 


onoftnnnnftf)  n^rioor* 


H2I I )=CAPAR( I ) /SV( I 1 
GO  TO  169 

167  M2< I1=CAPAC< I  l/SVt  1 ) 

169  HI  I )=CAPAR(  Il/SVII) 

H3I I )=CAPAR{ I l/SV( I 1*0£LTARI I ) 

031=G31+H3U1 

Q3I I+l )=Q31 

HU)=0.5*HU) 

H2( I)=0.5*H2(I> 

H31I1=0.5*H3II1 

ZERO  D1FF  INDICATORS  AND  X2 

X2 (I 1*0.0 
X3(I1=0.0 
X4II1=Q.0 
¥21 11=0.0 
TGI 1 1=0.0 
SUMX3I 11=0.0 
SUMX41 I )=Q.C 
PH1( I ) =0.0 
PHC1 11=0.0 
PR1 1 1=0.0 
152  RH0(I1=0.0 
X2UMP11*0.0 
PH01 IHP11=0.0 
PHI ( IMP1 1=0.0 
X31IMP11=0.0 
X4( IMP l 1=0.0 

X61IR+U=0.0 


X5(IR+11*0.0 

TG( IN1=0.0 
Y2< IN) =X6l INI 
X3 1  IN) =-1.0 


SET  UP  SOURCES 


*********************************************************************** 


FORM  Y2  ANO  TG 

SET  X3=— 1  IF  DIFFUSION  CRITERION  MET  USING  HCB 

00  156  I =1 N» ICX 
IF  (I-IM)  135,136,136 
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n  r>  n  nnr.nnn  on 


0 


k.;- 


135  T EHP 111 *H3 (I*1)*H3(I) 

IF  IAMIN1IH3II),H31I«-1))-CMIN)  U$7, 1157,1156 

1157  IFIABSI IH3I I 1-H3I I+1) ) /TEMPI  1) J— GA)  1156,1156,1158 

1158  IF  CX6U)  *EQ.  0.0)  60  TO  1159 
IFIA8S1IX6II)-XSII+1))/IX6II)+X61I+1)))>GA)  1156,1156,1159 

1159  TGII*1)*0.0 
GO  TO  156 

1156  TGIl4l)=lX6I I+1I—X6I I ) ) /TEMPI II 

TGI  I+1)*IQ1II  +  1)-QHI1)/IEMPI1I*IX511*1)+X5I  IU/2.0 

Y2I l6l)*IX6I I+1)PH3II ) +X6I I )  *H3(  1+1)  )/T£MPU) 

IFIABSI TGI 1*1) l-HCB*Y2I 1*1) )  158,158,156 

158  X3H6lJ*-i.O 
156  CONTINUE 

RADIATION  80UN0RY  CONDITIONS  AT  ICX 

136  Y21 JMP1)=X6I ICX) 

TGI IMP1)=0.0 

ENLARGE  TRANSPORT  REGION  80UNDRIES  BY  HVB 
MEAN  FREE  PATHS 

137  I-IN+1 

159  IF 1X311)1160,161, 162 

160  l»I+l 

IFI I -ICX-1 ) 159, 161 ,172 

162  SW3.162 
CALL  UNCLE 

161  J=I-1 

165  IFIQ3II)-Q3U)-HVB)  163,163,164 

163  X4IJ)=-l.O 
J*J~1 

IFI J-IN* 164, 165,165 

164  i*m 

IFlI-ICX-1) 166, 172,172 

166  IF(X3(I))168,164,162 
168  J*I 

171  IFIQ3I JJ-Q3I  I-D-HVB)  170,170,160 
170  X4I J) =— 1.0 

.)*J+1 

IFI J-ICX-l) 171,160, 160 

172  I*IN+1 

TEST  TO  FORM  TRANS  REGIONS 

IFIX3I 1N))173,174,162 

173  IFIX4I  IN)H74,175,162 

174  I AX=IN 

176  IF<X3(I))177,178,162 


67 


onno  non  oc>f>  or»o 


177  IF(X4(I1)176, 179,162 


REHOVE  ONE  ZONE  OIFF  REGION 

178  1«I+1 

Ifil-ICX-l) 176, 180,180 

179  i*i+i 

IF(I— ICX-1) 187, 180, 180 

187  IF (X3 (11)188, 176, 162 

188  IF( X4( I ) ) 176,189,162 

189  I8X=I~3 
GO  TO  181 

175  X2 ( IN)=C.O 

PR  ( IN) *Y2( I N)*45« 666 
RK0(IN)=Y2(IN)*137.0 
PH0(IN)=Y2<IN)*l,£i2 
PH1( IN)=PH0( IN) 

185  IF ! X3( I ) ) 182, 183, 162 

182  IF ( X4( I ) ) 183, 184,162 

FORM  X2  FOR  OIFF  ZONES  IN  ORDER 

184  X2( I  )--4. 12E12*TG< I )*CSQD( I ) 

PR  ( 1 I =Y2( I ) *45.666 
RH0( I ) =Y2( I ) *137.0 
PH0( I )=1*£12*Y2< I ) 

PH1( I)=PH0( I i 
1  =  1  +  1 

IF(l-ICX-l)  185,186,186 

00  TRANS  TO  IM  IN  REGION  OF  NO  SOURCE 

183  I AX= I 

GO  TO  178 

180  IBX=IM 

181  CALL  TRANS! IAX,IBX) 

IF(IBX-IM) 190,760, 760 

190  1*1 BX+2 
GO  TO  184 

186  IAX=I 

GO  TO  180 


* 

GPTIONAL  EDIT  OF  X2  ETC  * 

* 

C*********** ************************************************************ 
760  IF  (EDITMF)282,282,280 
280  IARG1=S0LI0( 181+0.001 
I ARG2=EDITMF+0«001 
IF  (M0D( IARG1 ,  IARG2) )  282,281,282 


r>ooo  o  o  o  finonn  oc*r>r>o 


281  WRITE  (3)  HNU*  IN,  IN,  IMP1,  SOLI0U8),  TH,  DHNU 

WRITE  (3)  (C(l)t  1»IN,  IMP1I ,  (H3(I),  I=IN,IHI,  (X6U),  I>IN,IMP1), 

1  IX2II),  I  *  IN, IMPI) ,  i PHO III,  I*IN, IMPI ) ,  (PHl(I),  I*IN, IMPI), 

2  (RHO(I),  I  =  IN,  IMPI) ,  ( PRI II,  I*IN»  IHP1),  (CSQD(I),  I*IN,!MP1) 
XX=-2.G 

WRITE  (3)  XX, XX, XX, XX, XX, XX, XX 

BACKSPACE  3 

JNUCT-1.0 

282  00  292  I=IN,IMP1 
SUMX2CI) =SUMX2(I )+X2< I > 

ERR ( I) »£RR ( I )+RHO{ I ) 

PRR ( I ) »PRR( I 1+PR< I ) 

292  CONTINUE 

ADVANCE  FREQUENCY  AND  STORE  EMERGENT  FLUX 
TEST  FOR  COMPLETION  OF  GROUPS 

HNUP*HNU 

HNUP4*HNU4 

CRTPC! IHNU+50)  *  X2(  IMPI  I 
IHNU  *  IHNU  ♦  I 

IF  (IHNU  -  NHNUi  283,  283,  285 

283  CALL  DVCHK  (KOOOFXI 

GO  TO ( 284, 701) ,KOOOFX 


ENO  FREQUENCY  LOOP 


* 

* 

♦ 


284  SI  »  13.0284 
CALL  UNCLE 

285  DO  286  1*INM1,IMP1 
288  X24 I )=SUMX2( I ) 

286  ER(  1 1  *  SUMX2U)  -  SUMX2U  +  1) 

FORM  MONOFREQUENCY  QUANTITIES  AND  FIND  MIN  TIME  STEP 

WSB»0.0 

DO  81  I=1,MAXLM 
81  WSB=WSB+ELM(  I ) 

DTR1=1.£10 

DTR2=1.E10 

FOR  MULTIFREQUENCY  WILL  CALL  KAPPA  ONCE  MORE  FOR 
AVERAGED  OPACITIES 

IF  ( KMAX  .£Q.  0)  GO  TO  1110 
CALL  KAPPA! IN, IM) 

1110  00  124  I=IN,IM 


c 

C  IF  ROSS  IS  ZERO  EXIT 

C 

IF  I CAP AR 1 1 ))illlfllllt1112 

1111  Sl=13.1111 
CALL  UNCLE 

C 

C  THICK  -  THIN  ETC 

C 

1112  TEMPI  3 )  =CAPARI 1 1 

IF I  SOLID! 10) ) 1113, 1114, 1113 

1113  TEMPI 1 I^CAPARI  1 1 
GO  TO  1117 

1114  TEMPI 1 )*SQRT (CAPARI I ) ♦CAPAC (III 

1117  IF I  TEMPI 11)1111,1111,1118 

1118  HI  I )=0»5*TEMP( 1 l/SVI I ) 

H3I I )  =  Hill  *  DEL TAR  1 1 ) 

IF  1 0.0  5-THETA  I  ID  1120 ,1201 , 1201 
1120  IF  I H3 1 1 1  .GT.  0.1)  GO  TO  1201 
IF  IERII)  .£Q.  0.0)  GO  TO  1201 
M$BB*E 1 1 )*GI 1 ) 

IF  ITELMI37)  .EQ.  0.0)  GO  TO  1201 
IF  (WSBB-TELMI37)*WSB)  1201,1203,1203 
1203  TEMPI 2)=0. 5*CV II)*THETAII)*GI I l/ABSI ERI I ) I 
GO  TO  1202 

1201  TEMPI2)  =  1.5  +  1.5  *  H3II)**2)  *  CVII)  /  I AC*TEMP(3)*THETAI I )**3) 
TEMP(2)=TEMPI2)*TELMI25) 

C*********************************************************************** 

c  * 

C  FINO  MINIMUM  TIME  STEP  * 

C  * 

C********>M***** I******************************************************* 

1202  IFITEMPI2II124, 124,1211 
1211  CONTINUE 

IF  I  TEMPI  2)-DTRi)U22, 1130, 1130 
1122  OTR2=OTR1 
IMN2=IMN1 
DTR1=TEMPI 2) 

IMN1=I 
GO  TC  124 

1130  IFI TEMPI  2) -0TR2) 1140,  124,124 
1140  0TR2=TEMP( 2 ) 

IMN2=I 

124  CONTINUE 

DTRM  IN=OTR 1 
£0= I MN 1 
C 

C  PRINT  MINIMUM  TIME  STEPS  BETWEEN  EOITS 

C 

IFIOTRl-TELM(26))jtl50,1160fll60 
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o  o  o  o  o  r>  r> 


1150  TELMI 26) =DTRi 
TELM(27)  =  IKM 
T£LM(28)=DTR2 
T£LMf29)=IMN2 
TELM!3C)=SGLID( 18)+ 1*0 
1160  CONTINUE 

DETERMINE  IF  RADIATION  OR  HYDRO  WILL  SUBCYCLE 

IF  { OTRMIN-DTR )  125,24G#1165 
1165  BLANK3*TH+AMIN1(0TRMIN,GR*DTH2) 

IF  I S17)  240  t J  23# 240 
123  S9  *  1.0 


* 

REDUCE  TIME  STEP  * 

♦ 

£*********************************************************«************♦ 
125  NRAD-ZP1 ( 18) /DTRMIN+1.0 
OTR*ZP 1 1 18) /FLOAT! NRAO ) 

IFINRAD-NTIMES)  240»240»238 
238  51*13.238 
CALL  UNCLE 
240  RETURN 
END 


f 


$I8FTC  STRANS  LI  ST, DECK, REF 
SUBROUTINE  TRANSIN, Ml 
C  OATE  COMPILED  3  JUNE  1965 
C  ANGULAR  INTEGRATIONS  ON  LINEAR  FORM 
C  CONNECTED  LINEAR-QUADRATIC  INTERPOLATION  AT  X=0. 
C  EDITING  OF  VARIOUS  M.F  - — 


v  cuiumi  ur  vakiuuo  n.r.  QUANTITIES  ON  S12  CYCLES 
C  TRDfiG  USED  TO  TRANSFER  FROM  THIN-PLANCK  TO  THICK-RQSSELAND 


C* 

c* 

c 


S  P 

UTTER  CO 

M  M  0 

N 

** 

♦ 

COMMON 

LMDa<37), 

NR  ,  NSMLR  , 

1A 

t 

IB 

,  1CA 

9 

1CB  , 

1 

KMAX 

,  8LANK1, 

BLANK2,  8LANK3, 

I  API 

t 

IBP1 

,  1CAP1  , 

1CBP1  , 

2 

II 

»  IG 

9 

NR AD  »  BLANK4, 

IAM1 

9 

IBM1 

,  ICAM1  , 

1CBM1  , 

3 

IIP1 

,  1GM’ 

9 

I ALPHA,  BLANKS, 

TH 

9 

TMAX 

,  BLANKS, 

DELPRT , 

4 

FREQ 

,  CNTMAX, 

AR  ,  ASMLR  , 

PUSHA 

t 

PUSHe  ,  B01LA  , 

B01LB  , 

5 

CVA 

,  CVf 

f 

SLUG  •  ALPHA  , 

HVA 

t 

HVB 

,  HCA 

9 

hcb  • 

6 

EMINA  ,  EMKN8  , 

CA  ,  C8  , 

GA 

9 

GB 

,  GL 

9 

GR  , 

7 

RHOL 

,  RHOR  , 

EPIO  ,  BPSl  , 

RIA 

t 

RIB 

,  RDIA  , 

RDIB  , 

8 

RPI A 

•  RPI8  , 

RPDIA  ,  RPDIB  , 

TPRINT, 

TA 

,  tb 

9 

TC 

COMMON 

TO 

9 

TE  ,  DTH2  , 

DTH2P 

t 

DTH1 

,  DTRM1N, 

DTMAX  , 

1 

DTMAX1,  DTMAX2, 

DTMAX3,  DTR  , 

SWITCH, 

CO 

,  CMIN  , 

DELTA  , 

2 

GAMA 

,  NCR  IT  , 

SIGMAQ,  AC  , 

AC03T4, 

CNVRT  ,  SUMRA  , 

SUMRB  , 

3 

ROI A 

,  ROIAM1, 

ROI 8  ,  ROI BP 1, 

GMS 

9 

SI 

,  S2 

, 

S3 

4 

S4 

,  S5 

, 

S6  ,  S7  , 

se 

9 

S9 

,  S10 

» 

Sll  , 

5 

S12 

,  S13 

, 

S14  ,  S15  , 

S16 

u 

S17 

,  S18 

, 

S19  , 

6 

S20 

t  EO 

, 

FO  ,  TAU  , 

ZERO 

f 

R 

(152), 

D£LTAR( 152) , 

7 

ASQ 

(1521, 

RD 

(152),  VD 

(152), 

ROD 

(152), 

SMLR 

(152), 

8 

OELR 

<  371, 

P 

(152),  PI 

(152), 

PB 

(152), 

PB1 

(1521 

COMMON 

P2 

( 152) ,  SV 

(152), 

RHO 

(152), 

THETA  (15 2), 

1 

M 

( 1 52 i , 

E 

(152),  El 

(1521, 

EK 

(152) , 

A 

1152), 

2 

V 

(152), 

G 

(152),  D 

(152), 

C 

(152), 

X2 

(152), 

3 

X3 

(1521, 

X4 

(152),  X5 

(152), 

X6 

(152)  , 

X7 

(152), 

4 

SMLA 

(1521, 

SMLB 

(152)  ,  SMLC 

(152), 

SMLD 

(152), 

SMLE 

(152), 

5 

EC 

(1521, 

ER 

(152),  SMLQ 

( 152)  » 

SMLH 

(152), 

BIGA 

( 152), 

6 

81 GB 

(152) , 

CV 

(152),  BC 

(152) , 

BR 

(152)  , 

CHIC 

(152), 

7 

CHIP. 

(152), 

CAP AC  (152),  CAPAR 

11521 , 

CRTC 

(152), 

CRTR 

(152), 

8 

CRTpc 

(152), 

GGFR  (152),  FEW 

(152), 

CAR 

(152) , 

OKLM 

(  37) 

COMMON 

TELM  (  37),  EKLM 

(  37), 

ELM 

(  37), 

FCLM 

(  37), 

1 

FRLM 

(  37), 

WLM 

(  37),  QLM 

(  37), 

AMASNGJ  37), 

CHRNO  (  37), 

2 

ZP1 

(  37), 

ZP2 

(  37),  SOLID 

(  37), 

£ 

CHCK 

(  37), 

RK 

(104), 

3 

RL 

<  37), 

RHOK 

(104),  RDK 

(104), 

THETAKI 104) , 

TEMP 

(  16), 

4 

HEAD 

(  12), 

MAXI 

,  MAXLM 

C*  ** 

Q***************  **************************#******.*******************:**** 

r  * 

COMMON  /UNDLY/  HNU,SGNL  ,IHNU,NHNU»HNUP,NT,  IM,  IN,  DHNU, THICK, NY 

COMMON  /CNTRL/  SCYCLE,  JMULT 


ws&ztgis  mm  was  mask  mm&m s  was  so?  pximsp  ' 
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f- 


k 

f 


l 


COMMON  /DAVIS/  FEX  (152)  ♦  JCX,  IR,  PR(152), 
1  TGI152) *  XC2400),  0X1152) 


DIMENSION 

CSQO 

11) , ERRlli 

f 

FH 

III, 

H 

ID, 

H2 

Hi 

1 

H3 

Hit 

H4 

ill,  PRR 

ID  * 

Q1 

ID, 

Q2 

(lit 

Q3 

<D 

2 

Q37 

m. 

Q38 

lilt  SUNX2 

UK 

SUKX3  111, 

SUHX4 

CD, 

X8 

ill 

3 

Y 

Cl). 

Y2 

(U,  YSQO 

m. 

PHO 

<11, 

PHI 

ID 

C 

EQUIVALENCE  C ACC3T4.TR0BG) ,  <BC  » SUHX4) e  (  BIGA  »Y 


1 

(3IGB  ,  H 

i  f 

(BR  ,K3  )» 

ICAR 

»Q37 

), 

(CHIC 

•SUHX3) , 

2 

CCHIR  ,Q38 

it 

ICRTR  , SUMX2) * 

(EC 

,Q1 

)• 

(GOFft 

,Q3 

i* 

3 

(ROD  *ERR 

it 

(SI2  ,EDITMFD 

( SMLA 

»H4 

J, 

(SMLB 

,  FM 

it 

4 

(SMLC  ,CSQD 

it 

( SMLD  ,Y2  i. 

CSKLE 

»H2 

), 

4SMLH 

,X8 

it 

5 

ISHLR  ,PRR 

K 

IV  ,02  i. 

(X7 

»YSQD 

it 

ICRTC 

tPHl 

it 

6 

CCRTPCrPHO 

) 

c 

* 

> ********************************** 

* 

C 

►♦A  *  |rr  t 

-r  ^ 

c 

ox 

CONTAINS  X  FROM  THE  PREVIOUS  Y  LINE 

* 

c 

♦ 

c 

CSQD 

SAME 

AS 

SKLC 

♦ 

c 

EDITMF 

SAME 

AS 

$12 

* 

c 

ERR 

SAME 

AS 

ROD 

* 

c 

FM 

SAME 

AS 

SMLB 

* 

c 

H 

SAME 

AS 

8IGB 

♦ 

c 

H2 

SAME 

AS 

SMLE 

* 

c 

H3 

SAKE 

AS 

BR 

* 

c 

H4 

SAME 

AS 

SMLA 

* 

c 

PHO 

SAME 

AS 

CRTPC 

* 

c 

PHI 

SAME 

AS 

CRTC 

* 

c 

PRR 

SAME 

AS 

SMLft 

* 

c 

Q1 

SAME 

AS 

EC 

* 

c 

Q2 

SAME 

AS 

V 

♦ 

c 

Q3 

SAME 

AS 

GOFR 

♦ 

c 

Q37 

SAME 

AS 

CAR 

* 

c 

Q38 

SAME 

AS 

CHIR 

* 

c 

SUMX2 

SAME 

AS 

CRTR 

* 

c 

SUMX3 

SAME 

AS 

CKIC 

* 

C 

SUMX4 

SAME 

AS 

BC 

* 

c 

TRDBG 

SAME 

AS 

AC03T4 

* 

c 

X8 

SAME 

AS 

SMLH 

* 

c 

Y 

SAME 

AS 

BIGA 

* 

c 

Y2 

SAME 

AS 

SMLD 

* 

c 

YSQO 

SAME 

AS 

X7 

* 

c 

* 

£**#**#*«*  $^**!*#****iJi*fc¥*****£:**$***££$***1tt**<'*****$I:<‘*****As(‘***:f***ifc***c* 

c 

* 

c 

PHO 

SIOE  FLUX 

(  £RGS/CM**2-5£C ) 

* 

c 

PHI 

FORWARD  FLUX 

i  ERGS /CM* *2- SEC) 

* 

c 

PH2 

BACKWARD 

FLUX 

(ERGS/CM**2-SEC) 

* 

f 

j 
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C  FRR  RADIATION  PRESSURE  I£RGS/CM**3I  * 

C  RHO  RADIATION  ENERGY  1ERGS/CM**31  * 

C  * 

q***^****?****^#*****#-?************#******************#*********  ****** 

c  * 

C  HCI)=G.5/ROSS£LAND  HFP  * 

C  TEMPI 1)=D£LTAX*H4I J  * 

C  H2C I )=0. 5/ PLANCK  MFP  * 

C  T£MH2I=DELTAX*H2U)  * 

t  H3< I )=i TAU  SUB  R)  * 

C  * 

£********-?*******«********************************#***************♦***** 

c 

IAX~N 
13  X=M 
SO  99 

99  Q2{ t )=0» 

CALL  OVCHK(KCOOFX) 

GO  T0a00,I01)tKG00FX 

100  51*14,0100 
CALL  UNCLE 

101  I8XP5«!BX-*1 
IALPHA=ALPHA 


ERROR  IF  NOT  SPHERE 


GO  TO  I 102,102,103) 
102  51*14.0102 
CALL  UNCLE 


I  ALPHA 


SPHERE  ONLY 

103  1=1 8XP1 
F2*C.O 

113  IFUBX-IM)  106,108,110 


ERROR  IF  INOEX  EXCEEDS  NORMAL  RANGE 


110  Sl=14.0ilO 
CALL  UNCLE 

£********#*****♦*#********#***$***##*******^<1**************  ************ 

c  * 

C  CALCULATE  Y  =  Q  RAY  * 

C  * 

C**************************************.#********.************************ 

c 

C  LHS  OF  RAY  FIRST,  STORE  F2  IN  SUMX3. 

C 

111  TEMP { 2 ) =H2I I ) *0£LTAR< I ) 

IF  I  TEMPI  2 )-fR08G)  131,131,132 


7.7 


non  non  nnn  non  non  nnn  nonnn  non 


SMALL  OPTICAL  DEPTH  EXPANSION 

131  F2*F2-MIY2ilHY2II  +  m*0.5+X6m-F2-F2**TEHP{2) 

GO  TO  112 

GENERAL  EXPRESSION  FOR  F2IY*0*) 

SAVE  EXP 

132  H4I I ) *FR£XPI— H3i Ill 
F2*YZIIHTGlIHCIF2-YZ(mi-TGCHl)> 

1  *H4£ I )+TGI 1*11— TGI  I ) }*H4I I J 

112  IFIF2.LT.0.1  GO  TO  115 

SAVE  LHS  INTENSITIES  IN  SUMX3 

114  SUMX3I I )=F2 

IFITGt  U.EQ.O.)  Y2U)*X6II-1) 

1*1-1 

IFITGC  I)  «£Q«0«  )  Y2II)=X6U) 

IF( I-IAX>  116«1I1*111 

F2-0 •  IF  1  GREATER  THAN  ICX 

105  1*1-1 

108  SUNX3I 11=0,0 

IFII-l-ICXJ  110,114,105 

DIFFUSION  BOUNOARY  CONDITION  AT  IBXP1 

106  IFITGI  D.EQ.O.)  Y2(I>=X6Cn 
TEMP(2)=H2I I  HDELTARI  I > 

IF  I TEMPI2 1-TRD8G)  1.33,133,134 

SMALL  OPTICAL  DEPTH  EXPANSIONIDIFF) 

133  F2*Y2I  I+D+TGI I  +  1H-II Y2II J  +  Y2I  HI)  H0.5+X6I  I)  )*T£HP(2) 
GO  TO  112 

GENERAL  EXPRESSION  FOR  F2IDIFF)Y=0. ,LHS 

134  H4I I ) =FREXP( — H3I I i ) 

F2  =  Y2II)+TGII)+ITGtI*l)-TGII)) *H4I I I 
GO  TC  112 

NEGATIVE  F2  ZEROED 

115  F2=0 .0 

GO  TO  114 
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uwu  ooo  UOU  O  U  O  U 'J  u  uuu  O  O  O  UUU 


RHS  OF  RAY 

116  1*IAX 

IAXP*IAX 

IF ( I AX-I NJ110 • 222 , 1 19 
222  IAXP*IN+1 

RHOi 1) *SUHX3! 11*137.0 
PR  (1J=SUNX3!1)*45.66 
PHC!1)*1.026E12*SUNX3!11 
PHII1I*PH0(1) 

60  TC  223 


0IFFU510N  BOUNDARY  CONDITION  AT  IAX 

119  TEMP(2I *H2! I-l)*DELTAR( I— II 

IF  <TEMPC2J-TRDBG)  135,135,136 

SMALL  OPTICAL  DEPTH  EXPANSION! DIFF) 

135  F2«V2U-ll-TG<I-l)*UY2UI*Y2II-lll»0.5  +X6!  1-1)  )*TEHPI2I 
GO  TC  121 


GENERAL  EXPRESSION  FOR  F2  CD1FF, Y«0*RHS i 

136  HA! 1—1 J=FREXP!-H3! I— 1) ) 

F2*Y2!  I  l-TG!  I  !♦!  TGC  II-TGU-l  I  l*H4I 1-1 ) 

GO  TO  121 


REGULAR  INTEGRATION  STEP  CY*0,RHS) 

120  TEMP<2I«H2!I-U*0ELTAR(I-1) 

IF  ! TEMP12I-TRD8GJ  I37,137f138 

SMALL  OPTICAL  DEPTH  EXPANSION  FOR  F2 

137  F2«F2*!(Y2UI*¥2<I-l)M0.5*X6!l-l)-F2-F2)*TEMPI2) 
GO  TO  121 


GENERAL  EXPRESSION  FOR  F2!Y*0«RHS) 

138  F2»Y2m-TGIIJ*<(F2-Y2!I-l)  +  TGU-m«H4U-l) 
1  ♦TG(n-TGU-lH*H4tI-l) 

121  IF(F2.LT.0.)G0  TO  221 

SAVE  RHS  INTENSITIES  IN  SUMX4 

223  SUMX4! I ) =F2 
OX! I ) =C! 1 1 

IFCTGCU.EQ.O.)  Y2!I)=X6m 
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-’•o 


1*1+1 

IFITGID.EG.O. )  Y2Ul=X6Ii-l) 

IF( I- IBXPI) 122,122,126 
122  IFU-l-ICX)  120,120,124 
C 

C  NEGATIVE  F2  ZEROED 

C 

221  F2*0-C 
GO  TO  223 

124  IFIF2.EQ.G.J  GO  TO  223 

F2=F2*FR£XPI-H3II-lJ-H3II-iJ I  *  * 

GO  TC  223 

C**»*»*»**»*»*«**»4*<*»*****»*»»»*»*»**»»»»«*»***»*»*»»»»»Me»*»*»»»»*»»» 
C  * 

C  COMPLETE  Y*0  RAY  * 

C  * 

126  YSQ01—0.0 
JJ=1 
JJJ=1 
KX*1 

C  * 

C  SEARCH  FOR  Y— LINE  NEAR  HALF  OF  C(IAX) 

c  ‘  •  1  ?  , 

C1=O.50*C<IAX) 

1I«IAX+1 

GO  TO  708  ’:**••• 

£***♦**♦************ ******************************************* ********* 

c 

C  SET  UP  Y  —  LINES 

C  * 

C*«*#****v** ************************************************************ 
C 

C  FIRST,  TEST  IF  Y  LIES  OUTSIDE  OF  ACTIVE  MESH. 

C 

127  IFIJJ— NY }  129,129,308 

128  IF(CCm-Y(JJ)I152, 150,150 
152  11=11+1 

GO  TO  128 

129  IF(C( ICX)-Y(JJ) 1307,151,151  i 

146  IFICUAX)-Y(Jjni28, 145,130 

130  IF(Ci-Y<Jjni4ia41,142 
142  IFICl-Yi  JJ+li  )  144, 140 » 140 

140  JJ*JJ+1 

KK *KK— I F I X f  X { KK+ 1 ) )+2 
GO  TO  127 
C 

C  SEARCH  FOR  Y-LINE  NEAR  THREE  QUARTERS  OF  CilAX) 

C 

141  C1=C1+C.25*C( I AX) 


*  * 


60  TO  145 

144  1F1V 4 JJ)*Y(JJ+1)-C1-C1) 140 ,141,141 

145  YSQDP=-XIKX| 

K*KK*1M-IBX*2' 

TEMP  45) =YSGDP- YSQDl 
60  TO  503 

151  IF4C4 I3XP1 )— Y4JJ I ) 300,300, i46 
C 

C  SEARCH  FOR  Y-LINE  NEAR  CUAX1 

C  .4. 

ISO  IF4C4 1  AXi-Y{  JJ.J-1) 1  153,153,145 
153  TEMPI  i)*X6|  II  H-xAUI?'!);- 
t  IF4 TEMPI ij.EQ.O.JGO  TO  154 


c 

c 

c 

c 


SKIP  Y-LI»iES  <UP.T0=4i  4f-^0URti  \x6l  AND  MFP  IH)  i££ 
FLAT  (CRITERION  BASED  ON  Gtl 


, IF4A8SI 4X6(11 1-X6U 1-I) l /TEMPI lj JrSt )  ,  l&Z  145*  145 
154  IF(AaS((H(IU^H(II-l)  )/(H(  li)*H(  H-l)  ))-eti  l55,145,145 
1§5  IF4JJ-JJJ-4J 157, 145,145 
157  j*l-  ,  .  , 

KK*XK-lFIX(X(KK+l))+2 
GO  TO  127 

307  TEMPI1 J=CSQD( ICX+1 ) 

GO-  ro  305 

308  TEMPI 1)=CSQD( IR+ll 

.-*•  GO  TO  305  -  . .  ^  _ _ 

£*#*** ********** ♦♦************#*************»***»*********i*******A***i* 

c  ......  .c,  ;* 

C  COMPLETE  X2  INTEGRATION  WHEN  LAST  Y-LINE  USED  * 

C  •  .  w  _ _  .1 

C^***********#***************************^*****^***********#********** 
;  300  DO  310  I=IAXP,IBXP1 


c 

C 

c 

c 

c 

c 


ANGULAR  INTEGRATION  OF  X2  USING  LINEAR  INTERPOLATION 
X2  ( I )  =X2  ( I)  +  (0X1 1 )  4-OX  { I  ))  *  {  SUMX’4{  I  )  -  SUMX3II I )  *0X41 J 
OPTIONAL  CALCULATION  ON  Si‘2 


IF4EDITMF)310, 310,625 
625  IARG1=S0LID(18)+0.001 
IARG2=EDITMF+0.001 
IF (MODI  I ARG1, 1 ARG2) >310,627,310 
627  IFU8X-IM)  215,317,510 
315  IF( I-IAXP )  510,312,340 
312  XS=SQRT(CSQD(  IB.XP1+1)-CSQD(  IBXP1M 

FMUBXPl)=Y2U8XPli+XS/C(IBXPl  +  l)*TG(  I8XP1+1)* 
l  EXP{-XS*H{IBXPl)i 


GO  TC  340 
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317  FMIIBXPU*0.O 
340  TEMPI51=CSO0UBXP11-YSQO1 
TEMPI 111*  CSQO I I J-VSQDi 

FU* I TEMPI  11 1*1  FMI 1  BXPll *FMII  8XP11 1  *1 CSQOI I8XP1I- 
1  CSQOI I) 1*1 SUMX4I IAXP)*SUMX31IAXP) I 1/TEKPI51 

FL*SUMX3I I 1+SUHX4I 1 1 
TEMP 161=0X1 I 1*0X11 ) 

FP=FL*FL 

?PL*FP*FL 

RHOII)  =P.HO(  1 )  *0X1  U*FP 

PR  1 I 1*PR| 1 1 *T£MPI61 *1 FPL *FL ) *0X1 1 1 

PHO (1 1 =PHO 1 1 1 *FU*( 0X1 II *TEMP 1 9 1 *CSQO 1 1 1 *ARS INI 0X1 1 ) /C 1 1 >) I 
1  *.667*IFL-FU)*1CII)*CSQD( 1 1-TEMPI 9)*YSQ0U  /0X1II 
PH11IJ=PHHI1*FPL*T£HP<61 
310  CONTINUE 

C******************************* ***♦*♦*♦♦*»**♦♦*♦*♦*♦»*♦*♦»*♦*♦*♦»*♦♦♦♦* 

c  * 

C  COMPLETION  OF  X  2  INTEGRATION  AT* 

C  ENO  OF  TRANS  REGION  * 

C  * 

C******************* **************************************************** 

c 

C  BRIGHTNESS  PRINT 

C 

320  IFICBJ  726,725,726 

726  IARG1=S0LI0( 18 1*0.001 

IF  IIARG1  .EQ.  01  GO  TO  725 
IARG2=CB+0.001 

IF  IM0DCIARG1,IARG21)  725,727,725 

727  XX  =  -1.0 

WRITE  (3)  XX,  IAX,  IBX,  NY,  S0LI0I18I,  TH,  XX 
WRITE  <31  I Y{ J) *  J=1,NY),  IQ2IJ.*,  J=l,NYl 
XX  =  -2.0 

WRITE  (3)  XX,  XX,  XX,  XX,  XX,  XX,  XX 
BACKSPACE  3 
725  00  321  1=1 AX, IBXP1 
C 

C  OPTIONAL  CALCULATION  ON  S12 

C 

IFIEOITMFl  321,321,730 

730  IARG1=SOLIO< 181*0.001 
IARG2=E0ITMF*0.001 

IF <MQO( I ARGi, I ARG2) 1321,731,321 

731  IF(I.EQ.l)  GO  TO  321 
RHOII  )=RHOI  1 1*34. 25/CU) 

PRI I1=PRII )*5.71/(CSQO( I)*CI I) 1 
PH01 1 1 =PHOI I )*0.327E12/CSC01 1 1 
PHI I I1=PH1C 11*0. 343E12/CSQ0I II 

321  X21 1 1=X2( 1 1*1. 026E 12 

322  RETURN 
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c 

Z  COMPLETE  X2  INTEGRATION  WHEN  NO  MORE  SOURCES 

C 

305  F2S0*0 

TEMP<51*TEMPI1)-YSQD1 
00  306  1=1 AXP, I8XP1 
TEMP£ 151=T£MPI 1 )~CSQO( I 1 
IFITEMPI 151)  318*319,319 

318  FMI I )=0.0 
C 

C  ANGULAR  INTEGRATION  OF  X2  USING  LINEAR  INTERPOUTION 

C 

XM*SQRT I— TEMPI 15)1 

X2I I 1-X2I l 1 *<0X1 1 l-XM) *1 OXI I 1 +0X1 1 )+XM)*ISUNX4{ I 1-SUMX3I 1 1 1 
IFIEDITMF)  306,306,323 

323  IARGl=S0LIDI18)+0.001 
IARG2=£DITMF«0.001 

IF (MODI IARGl , IARG21 I  306,324,306 

324  FL=SUMX3I I 1+SUHX4I I 1 
TEMPI 6 1=0X1 I 1*0X1 I 1 
TEMPI7)=-TEMPI 151 
FCA=SQRTI TEMPI  1)1 
FXM=0XU1-XM 
FXP=GXI I 1+XM 

RHO 1 1 1 =RHO 1 1 1 +FXM*FL 
PRI I )=PRI I )-FXP*l TEMP  1 61+TEHPI 7> )*FL 
1  +4.C*F.,„*GXI  1 1*TEMPI6J 

PHOI I 1 =PHO I i 1-IFL*XM*I OXI I 1 *TEMPI 9 )-XM*FCA+CSQDi I 1 * 

1  (ARSINIOXIll/CIIl  1-ARSIN1 XM/CI 1 1 1 1-'0.667*FL*I  FCA*TEMPI  1) 

2  -TEMPI 9)*YSQD111)/FXM 

PHI 1 1 ) =PH1 1 1 } +FXM*FL*I FXP+OX III) 

GO  TO  306 
C 

C  ANGULAR  INTEGRATION  OF  X2  USING  LINEAR  INTERPOLATION 

C 

319  X2U  >=X2 (11  +  ( OXIll+OXI I )  1*1  SUMX4I I1-SUMX3I  I)  )*OXI  I) 

c 

C  OPTIONAL  CALCULATION  ON  S12 

C 

IF! EOITMF) 306,306,635 
635  IARG1=S0LIDI 181+0.001 
IARG2-E0ITMF+0.001 
I F I MOD  1 1 ARG 1 , IASG2 1 ) 306, 63  7, 306 
637  TEMPI5)=TEMPI 1J-YSQD1 

FU=TEMPI 15)*ISUMX4I IAXPJ+SUMX3IIAXP) 1 /TEMPI  5) 

FL»SUMX3H  1+SUMX4I 1 1 
TEMP (6)=0X1 I ) *0X1 1 1 
FP=FL+FU 
FPL=FP+FPL 

RH04 I ) =RHO( I J+OXI I)*FP 


! 

1 


83 


nnr>  n  n  n  nnn  non  non  non  o  o  o  o  o  o  o 


PR< I  )=PRII  )*T£MP(6)*1FPL+FL)*0XU ) 

PH0< I I  *PMOII J+FU*( GXI I ) *T£MP(  9) +CSQOC II*ARSIN(OX( I )/C( I) i i 
1  +  .667*<FL-FU)MC(n*CSQQ(I)-TEMP(9)*YSQDl)  /  OXU) 

PHU  I )  *PH1  <  1 1  +FPL*T£MP  (61 
306  CONTINUE 
GO  TO  320 

£*♦*♦♦*♦*«*♦**♦♦#♦*♦ ************************************ **************** 

♦ 


TYPICAL  Y-LINE  INTEGRATION 


* 

* 


*********************************************************************** 


LHS  CALCULATION  FIR5T,ST0R£  F2  IN  FH 


503  1=1 8XP1 
F2=0.0 

502  IF( IBX-IMI  506*508,510 

ERROR  IF  INOEX  EXCEEDS  NORMAL  RANGE 


510  Sl=14.0510 
CALL  UNCLE 


DIFFUSION  BOUNDARY  CONDITION  AT  IBXP1 

506  IF  ( TGI  I )  .EG.  0.0)  Y2U)*X6(II 
X8 (11 =X(K)/C (I) *TG(I) 

TEMP(2)=(XtK— II— XIK))*H2(I ) 

IF  < TEMPI 2)-TR0BG)  533,533,534 

SMALL  OPTICAL  OEPTH  EXPANSION  (DIFF) 

533  F2=Y2(  l*il+X(K-U/C(I  +  l)*TG(  I+1J+4  (Y2(l)+Y2(  1+1)  )*0.5+X6( II) 
1  *TEMP(2) 

GO  TO  512 


GENERAL  EXPRESSION  FOR  F2  (DIFF, LHS) 

534  TEMP(i)=(X(K-l)-X(K))*H(I) 

H4( I )=FREXP(— TEMPI  I)) 

F2=Y2(  l  )+X8  ( I  )♦{  X(  K— 1 )  /CC  I-o-l  )*TG(  I+ll— X8<  I  ) J*H4(  I ) 
GO  TO  512 


SAVE  X8  FOR  RHS=(DIFF  INTENSITY) 

511  X9(I)=X(K)/C(I)*TG(I) 

TEMP (2)=(X(K-1I—X(KI)*H 2(1) 

REGULAR  INTEGRATION  STEP(LHS) 
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I 


f 

t 


IF  ( TEMPI2I-TRDBG)  531,531,532 
C 

C  SHALL  OPTICAL  DEPTH  EXPANSION  (LHS1 

C 

531  F2=F2+I I Y2(I 1+Y2I I+i) )*0.5+X61 I I-F2-F2)  *TEMPI2) 

GO  TO  512 

C 

C  GENERAL  EXPRESSION  FOR  F2ILHSI 

C 

C  SAVE  EXP  FOR  RHS 

C 

532  TEMPI1)=IXIK-1)-XIK))*H1I) 

H4I I )=FREXPI-TEMPI 111 

F2=Y2II)+X8II)  +  IlF2-Y2II*l)-X8U+l))*H4tI) 

1  +X8I I  +  D— X8I I))*H4III 

512  IFIF2.LT.O.)  GO  TO  515 
C 

C  SAVE  F3  OF  LHS  IN  FM  FOR  INTEGRATION. 

C 

613  FMII)=F2 

IF(TG(  D.EQ.O.)  Y2I I )=X6I I-i ) 

1=1-1 

K=K+1 

IFITGID.EQ.O.)  Y2I  l )=X6I I ) 

IF (I-I AX >516 ,513, 513 

513  IFIXlK))514,5C9,5ll 

514  TEMPI  31 =SQRT (DEL TAR ( I )*I4.0*CI I )+DELTARI I ) ) ) *HJ I ) 

IF  ( TEMPI3)  .GT.  0.0)  GO  TO  517 

S1=14.0517 
CALL  UNCLE 

517  TEMP 1 2)=XI K~l )*H2I I ) 

IF  ( TEMPI 2J-TR08G)  535,535,536 

535  TEMPI16)  =  IY2<I  +  i)+X6lII)*T£MP(2i+I0.667*Y2tI)-0.5*Y2IH-ll 

1  -0. 167*X6I I ) )*TEMPI3 )*H21I ) /HI I ) 

F2=F2*I 1.0-2. 04 TEMPI  2) )+T£MPI 16) 

GO  TO  521 

536  TEMPI 1)=2.0*X!K“1)*HI I) 

H4I I  I'-FREXPI-TEMPC  1)  ) 

TEMPI 7 ) =FREXP I— TEMPI  3) ) 

TEMPI13)=IX6II)-Y2II) ) /TEMPI  3) **2*2.0 
TEMPI  15)  =  1X61  D-Y2II  +  1)  ) /( TEMPI  1) -TEMPI  3)  ) 

F2=Y2< I)+TEMPI13)+TEMPI 7) *(— TEMPI 15)-TEMP( 13 )*I TEMPI 3)+l.C) ) 
1  +H41 1 1*1 F2- Y2 1 1* 1) ♦TEMPI 15) ? 

GO  TO  521 
C 

C  NEGATIVE  F2  ZEROED 

C 

515  F2=0.0 

GO  TO  613 
C 
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C  FIRST  TRANSPORT  ZONE  ON  LHS 

C 

5Q7  X8(n=X(K)/Ct  I)*TG(Ii 
GO  TO  613 

X  «  ZERO  ERROR 

509  $1=14.0509 
CALL  UNCLE 
508  FHII)=G.O 
1  =  1-1 
K*K+1 

IF{ I-l-ICX ) 510  » 507*  508 
523  F2=0. 0 
GO  TO  522 


DO  y-INTEGRATION  TO  X=0 

521  IFIF2.LT.0. )GO  TO  523 

522  FM< 1 )=F2 


FGRM  TOP  SLICE  CONTRIBS  TO  X2 

DO  524  J  =  I AXP  *  I 
333  TEMPUl)=CSQDiJ)-YSUDl 

X2( JI=X2(JI*(OXI JJ+UXI J) )*( SUMX4C JJ-SUMX3I J) )*GXT  J) 

OPTIONAL  CALCULATION  ON  S12 

IF(EBITMF) 524*524*331 

331  IARG1=SOLIDU8)+0.OJ1 
IARG2=EDITMF+0.001 

IF I NOD ( I ARG1 , 1 ARG2 ) ) 5  24 , 332 , 524 

332  FU»IT£NP(ll)*IFHII)+F2)+<YStiDP-CSQD<  J)I*(SUMX4I IAXP) 

1  +SUMX34 IAXP)))/TEMP(5) 

FL=SUMX3(J)+SUMX4(J) 

TEMPI6I=QX( J)*OX{ J ) 

FP=FL+FU 

FPL=FP+FL 

RH0( J J  =RHOI JI+OXI J ) *FP 

PRI J)=PR< J)*TEMP(6)*(FPL+FL)*0X(J) 

PH0( J)=PH0< J)+FU*<0X< J)*TEMPi 9I+CSQD4 JI*ARS1NI0X( J)/C< Jl)  I 
1  *-0. 667*(FL-FU)  *(G  (  l )*CSQDI I )~TEMP( 9J*YSQD1 >/OXt I  I 

PH1( Jl =PH1< J I +FPL*TEMP(6) 

524  CONTINUE 
I AXP=I + 1 
SUMX3I I J  =F2 
SUMX4I I )=F2 
1=1 AXP 
K=K- 1 
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IF  I TEMPI2I-TRDBG)  533*538,700 


SMALL  OPTICAL  DEPTH  EXPANSION  <X-0,RHS) 

538  F2*F2*(1.0-2*O*TEMP<2I)*TENP<16I 
GO  TG  704 

700  F2=Y2I  1 J+TEMPI  15HFREXPI-TEMPI  1 1-#-T£MP!3)  I* C— TEMPI  151 

1  +TEMPI 131*1 1.0-TEMPI 3) I )*H4I I“1I *IF2-Y2I I~l)— TEMPI 13) J 

GO  TO  704 


CALCULATE  RHS  Y-LINE  INTEGRATION 


516  IAXP=I AX 
I  =  IAX 
K=K- 1 

IFIXIK+1J ‘525,509,526 
526  X8U-l)=X(K+ll/CII-l)*TGU-ll 
TEMP<2)=IX(K)-XIK+1J )*H2«I-1I 

DIFFUSION  BOUNDARY  CONDITION  AT  IAX 

IF  I TEMPI2I-TRDBG)  537,537,702 

SMALL  OPTICAL  OEPTH  EXPANSION  IDIFF.LHS) 

537  F2*Y2<I~I)-X8(I-l)  +  {<Y2(n  +  Y2U-l)  1*0-5  ♦  X6I  I-1))*TEMP{ 21 
GO  TO  704 


GENERAL  EXPRESSION  FOR  F2IDIFF,RHS) 

702  TEMP (1I={X(K)-X{K+1))*H(I“1I 

701  F2=Y2IU-X8im(X8(I)-X8U-l)}*fREXP{-TEMP{l) ) 

GO  TO  704 

710  TEMPI2)=<X(KJ-X(X*1) ) *H2C I— 1 > 

IF  4 TEMP ( 2 i-TRDBG)  539,539,540 

SMALL  OPTICAL  DEPTH  EXPANSIONIRHS » 

539  F2=F2+UY2II  )  +  Y2U-l)  )*0.5+X6(  I-D-f  2-F2)*TEMP< 2# 
GO  TO  704 

540  F  2=Y2(  I )-X8II i<  < (F2-Y2I I~l) +X8( I-l ) )  *H4<  I-  11 
1  ♦Xfi(I)-X8(I-m*H4II-l) 

704  IFIF2-LT.0.IG0  TO  705 

FORM  CONTRIB  TO  X2 


706  CONTINUE 

ANGULAR  INTEGRATION  OF  X2  USING  LINEAR  INTERPOLATION 


8-7 


n  o  n  no  o  non  non 


703  X2III=X2<I1  +  I  (XlKI+X|Kl+GXIIJl*IF2~Fmill 

1  +  lGX(Il+0X(tj+X(K.J  1*ISUMX4I 11— 5UMX3IIJ  11*10X1  D— XIK1  J 

OPTIONAL  CALCULATIONS  OF  H.F.  QUANTITIES  ONLY  ON  EOITHF 

IFI ED I TKF 1735*735*736 
736  IARGl=SGLIDI181+0.00i 
I ARG2-ED IT MF+O .00 1 
IF  (NODI IARGIt I ARG211  735,732,735 
732  FU=FH<ll+F2 

FL=SUMX3t I J+SUMX4I I 1 
TEKP(61=0X(I )*QXII 1 
TEMPI7)=XI K1*XIK1 
FXH=GXI11-XIK1 
FXP=OXIIJ+X(K) 

FMI=Fli— FL 

RHO( I 1=RH0I ll+FXM*{FU+FLl 
PR ( 1 1 =PR (I I +FXP* ( T EMP ( 61 +T EMP I  7 1 1 *FH 1 
1  +4.Q*IFL*OXIll*T£HP{Gl-FU*XItO*TEMPC71) 

PHC(Il=PHCUl  +  IIGX{Il*FU-X{K|*FIJ*(GXUl*TEH!>{91-XtKl*Yl  JJ) 

1  +CSQDI I t*( ARSINIOXI It/Cil) 1-AFSINI XIR1/CI II ) 1 }-0.667*FHl 

2  *  I Y I  J«1 1  * YS  QOP— TEMP  4  91  *Y  SQD 111 /FXM 

PH1( II =PH1( I 1+FXH* (F2*IFXr+X(KJ J+S0KX41 I 1*CFXP+QX{ II 11 

SAVE  F2  AND  FM  FOR  NEXT  Y-LINE 

735  SUMX4I I )=F2 

SUMX3I I 1=FMI t 1 
GX(I1=X(K) 

IFITGU1.EQ.C.)  Y21 1 1=X6I  II 
1=1  +  1 
l 

IFITGI Il.EQ.O.)  Y2( 1 1 =X6II-11 
IF( I— I8XP1 1 707 , 707 ,706 
707  IF ( I -1-1 CXI  710,710, 711 

NEGATIVE  F2  ZEP.OED 

705  F2=0.0 

GO  TO  706 

711  IFIF2.EQ.0.1  GO  TO  706 

TEMP (ll=IX(K)-X4K+ll)*H(I— 1) 

H4( 1-1 1 =FREXP( -TEMP ( 1 1 -TEMPI  11 1 
F2=F2*H4( 1-11 
GQ  TO  706 

DIFFUSION  BOUNDARY  CONDITION  WHEN  X=0. 

525  X8I 1-1 1 =0.0 

TEMP ( 1 1 =X( K 1 *H( 1-11 
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50  10  7£l 
708  SHNU=H*iUP-HfcU 
712  IF<JJ.EQ»U  GG  TO  722 

£»»»»»♦»:»♦»»♦*»♦*.»  ♦♦***y**»  ♦»♦*»»*»***♦**»»**#  ■»**»*»**^*»*y*»»**»*»*»*»» 

* 

tiJiT  A  N  0  Y  ‘  ^  t  «  g  A3VA«Cc  * 

* 

C  ♦♦♦♦♦♦♦♦♦♦♦♦»  »***»*^»»*.*»»g»*******:fc»*5:g*****»*  *»*♦»»»****♦  »♦♦*»»»♦♦*♦* 
ftK-KK—  JFIXi  X IKK* II  )  *-2 

722  YS0C1=YSQ0UJJ 
TEHPm=Y[Jj} 

724 

Q2 1 J J~ i ) =f 2*3 . 274* 4 
SC  TO  12? 

END 
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Appendix  C 
FREXP 


fmnmi  Bias  «as  taesaOT  «as  w»  mts© 


U<,J4HWVXM*R»**V1 


${ 8HAP 
♦ 


* 

♦ 

♦ 

* 

* 

* 

♦ 

* 


FREEKA  LIST  »0£CK«R£F 

BECK  FREEMA  EVALUATES  FACTION  FREXP (ARG) 

EVALUATE  THE  NEGATIVE  POKER  OF  E<E**-X)FOR  X  BETWEEN 
9.0  AND  -10.0  TO  FOUR  SIGNIFICANT  FIGURES.  WHEN  X  LIES 
OUTSIDE  THE  RANGE,  THE  STANDARD  ROUTINE (EXP)  IS  USED. 

THE  HETHO0  IS  AS  FOLLOWS*  E**-X=£**-N*E**-P*E**-2 « 

WHERE  N  IS  THE  INTEGRAL  PART  Or  THE  EXPONENT  AND  P  IS  THE 
MULTIPLE  OF  2**-7  CONTAINED  IN  THE  FRACTIONAL  PART. 

R-X-N-P.  THE  VALUES  OF  £**-N  AND  E**-P  ARE  FOUND  BY 
TABLE  LCOX-UP  IN  ITA3  AND  FTAB  RESPECTIVELY,  E**-R  IS  FOUND 
AS  1.0— R. 


FREXP 


ENTRY 

EXTERN 

A  XT 

AXT 

TRA 

SXA 

LAC 

CL  A* 

CAS 

TNI 

TNI 

~TA 


UFA 

CAS 

TRA 

TRA 

PAC 

STQ 

SXA 

CLA 

UFA 

PAC 


FXP 

FAD 

STO 


FREXP 

EXP 


**,2 


RESTORE 

INDICES  2  AND  A 
RETURN  ADDRESS 
SAVE  INDICES  2  AND  4 
LOCATE  ARGUMENT 
PICK  Ur  THE  ARGUMENT 
SELECT  THE 
NUMBERS  BETWEEN 
0.0  AND  -9.99 
OUTSIDE  THAT  RANGE 
ARGUMENT  IN  THE  RIGHT  RANGE 
=0633000000000  EXTRICATE  THE  INTEGER  PART 
=0633000000000 
NORMAL 


** 

*-3,4 
FREXP ,4 
2,4 

=—9.99 

*+3 

*♦2 

NORMAL 


SIMPLE 
0,4 
T 1 

FREXP-2,2 
T 1 

=0624000000000 


PUf  THE  INTEGER  PART  IN  THE  INDEX 
PUT  THE  FRACTIONAL  PART 
SAVE  INDEX  2 
IN  THE  AC 

EXTRICATE  THE  2**~7  MULTIPLE 


0,2 


FTAB, 2 
FTAB, 2 
T 1 


PUT  IT  IN  INDEX  2 
REMAINDER  LEFT  IN  MO 
FIND  THE  FRACTIONAL 
POWER 
PUT  IT 


LOQ 

II 

IN  THE  AC 

FHP 

I TAS ,4 

MULTIPLY  IT  BY  THE  INTEGRAL  POKER 

TRA 

FREXP-2 

RETURN  WITH  ANSWER  IN  THE  AC 

• 

THE  EXPONENT 

LIES  BETWEEN  0.0  AND  -1.0 

SIMPLE 

STG 

11 

CLA 

T 1 

UFA 

=0624000000000 

PAC 

0,4 

FHP 

FTAB, 4 

FAD 

FTAB, 4 

- 

* 

TRA 

FREXP- 1 

vizcirms  east  ms 

m&m  msmtTm  ms  so?  mtm> 

NORMAL 

CLA 

2,4 

STA 

Ml 

TSL 

EXP 

TXI 

eci 

l.FREXP 

Ml 

pze 

♦* 

TRA 

FREXP-1 

T1 

PZE 

I  TAB 

OEC 

1.0 

0 

DEC 

0.3678794411 

1 

OEC 

0.1353352832 

2 

OEC 

0.0497870683 

3 

OEC 

0.0183156388 

4 

OEC 

0.0067379469 

5 

OEC 

0.0024787521 

6 

OEC 

0.0009118819 

7 

OEC 

0.0003354626 

8 

OEC 

0.0001234098- 

9 

FTAB 

OEC 

1.0 

0 

DEC 

.99221794 

DEC 

•98449644 

DEC 

.97683503 

DEC 

.96923324 

DEC 

•96169060 

OEC 

.95420668 

DEC 

.94678098 

DEC 

.93941307 

DEC 

.93210250 

DEC 

•92484882 

DEC 

•91765159 

DEC 

•91051038 

DEC 

,90342472 

DEC 

.89639421 

DEC 

.88941842 

DEC 

.88249691 

OEC 

•87562926 

OEC 

•86881506 

OEC 

.86205389 

OEC 

.85534534 

OEC 

,84868899 

OEC 

.84208444 

DEC 

.83553128 

DEC 

.82902913 

DEC 

.82257757 

DEC 

.01617622 

DEC 

.80982469 

DEC 

.80352259 

OEC 

.79726952 

OEC 

.79106511 

OEC 

.78490900 

94 


DEC 

•77880079 

DEC 

.77274011 

DEC 

.76672660 

DEC 

.76075939 

DEC 

.75483961 

DEC 

.74896540 

DEC 

.74313691 

DEC 

.73735377 

DEC 

.73161563 

DEC 

.72592215 

DEC 

.72027299 

DEC 

.71466778 

DEC 

.70910619 

DEC 

.70358787 

DEC 

•69811252 

DEC 

.69267976 

DEC 

•68728926 

DEC 

.68194076 

DEC 

.67663385 

DEC 

.67136824 

OEC 

•66614361 

DEC 

.66095965 

DEC 

.65581602 

OEC 

.65071242 

DEC 

.64564854 

DEC 

.64062406 

OEC 

.63563868 

OEC 

.63069209 

OEC 

.62578401 

DEC 

.62091413 

DEC 

.61608213 

DEC 

.61128774 

OEC 

.60653067 

OEC 

•60181060 

DEC 

.59712727 

DEC 

.59248041 

OEC 

.58786967 

OEC 

.58329484 

OEC 

.57875560 

DEC 

.57425170 

DEC 

.56978284 

DEC 

.56534874 

OEC 

.56094916 

DEC 

.55658382 

OEC 

.55225246 

OEC 

.54795480 

OEC 

.54369058 

OEC 

.53945954 

DEC 

.53526144 

OEC 

.53109600 

OEC 

•52696298 

DEC 

.52286212 

OEC 

.51879318 

OEC 

.51475589 

OEC 

.51075003 

OEC 

.50677535 

OEC 

.50283159 

OEC 

.49891852 

OEC 

.49503590 

OEC 

.49118350 

OEC 

.48736108 

OEC 

.48356841 

OEC 

.47980525 

OEC 

.47607137 

OEC 

.47236656 

OEC 

.46869057 

OEC 

.46504319 

OEC 

.46142419 

OEC 

.45783336 

OEC 

.45427047 

OEC 

.45073532 

DEC 

.44722766 

OEC 

.44374731 

OEC 

.44029404 

OEC 

.43686765 

OEC 

.43346792 

OEC 

.43009464 

DEC 

.42674762 

OEC 

.42342664 

OEC 

.42013151 

OEC 

.41686203 

OEC 

.41361798 

DEC 

.41039918 

OEC 

.40720543 

OEC 

•40403653 

OEC 

.40089229 

OEC 

.39777252 

OEC 

.39467703 

OEC 

.39160563 

OEC 

.38855813 

OEC 

.38553435 

OEC 

.38253409 

OEC 

.37955719 

OEC 

.37660345 

t)EC 

.37367270 

DEC 

.37076476 

OEC 

.36787944 

END 

SIBFTC  KAPPA  FULIST,0£CK,REF 
SUBROUTINE  KAPPA (N»M) 

C 

C  THIS  ROUTINE  NAS  RECOMPILED  JANUARY  II 9  1965 

C  TWG-OIANE-TAPE  YERSION 

C  ASSUMES  EXTRA  CALL  FOR  GREY  DATA  —  NO  GOOD  WITH  STANDARD  RAOTN 
C  DOES  MGNGFREQUENCY  WITH  AVERAGE  KAPPAS 

C 

c*********************************************************************** 


c 

‘  *. 

c* 

S  P 

U  T  T  E 

R  C 

0  M  M  0 

N 

- 

c* 

** 

c* 

.  ** 

COMMON 

LMGAI37J, 

NR.  , 

NSMLR 

♦  IA 

,  IB 

,  ICA  , 

ICB  , 

I 

KMAX 

,  BLANKlf 

BLANK2, 

BLANK3 

,  I  API 

,  I8P1 

,  ICAPi  , 

ICBPl  , 

2 

II 

,  IG 

t 

NR  AD  , 

BLANK4 

*  I  AMI 

*  IBM1 

,  ICAM1  , 

ICBM1  , 

3 

IIPI 

,  1GMI  , 

IALPHA, 

8LANK5 

»  TH 

,  TMAX 

,  BLANK6, 

OELPRT, 

4 

FREQ 

,  CNTHAX, 

AR  , 

ASMLR 

,  PUSHA 

,  PUSHB  ,  BOILA  , 

BOILB  , 

5 

CVA 

,  CVB 

9 

SLUG  , 

ALPHA 

,  HVA 

,  HVB 

»  HCA  , 

HCB  , 

6 

EMINA  ,  EMIN8  , 

CA  , 

CB 

,  GA 

*  GB 

•  GL  , 

GR  , 

7 

RHOL 

,  RHOR  , 

EPIC  , 

EPSI 

,  RIA 

t  RIB 

,  ROIA  , 

RD:?B  , 

8 

RPI A 

,  RPIB  • 

RPOIA  , 

RPDIB 

•  TPRINT,  TA 

t  TB  , 

TC 

COMMON 

TD 

9 

TE  , 

0TH2 

,  DTH2P 

,  DTH1 

,  DTRMIN, 

OTMAX  , 

1 

DTMAX1,  DTMAX2, 

OTMAX3, 

OTR 

,  SWITCH,  CO 

,  CHIN  , 

DELTA  , 

2 

GAMA 

,  WCRIT  , 

SIGMAQ, 

AC 

•  AC03T4,  CNVRT  ,  SUMRA  , 

SUMRB  , 

3 

ROI A 

,  ROI AMI, 

ROIB  , 

R0IBP1 

,  GMS 

,  SI 

,  S2 

S3  , 

4 

S4 

,  S5 

t 

S6  , 

S7 

•  S8 

.  S9 

,  S10  * 

Sll  , 

5 

S12 

,  SI3 

t 

S14  , 

S15 

,  S16 

,  S17 

,  S18  , 

S19  , 

6 

S20 

*  EO 

• 

FO  , 

TAU 

,  ZERO 

,  R 

(152) ,  DEL TAR( 1521 , 

7 

ASG 

(152), 

RD 

( 152) 

,  VD 

(152)  , 

ROD 

(152),  SMLR 

(152), 

8 

DELR 

1  37), 

P 

(152) 

,  PI 

(152), 

PB 

(152),  PB1 

(1521 

COMMON 

P2 

U52) 

t  sv 

(152), 

RHO 

(152),  THETA  (152), 

1 

W 

(152), 

E 

(152) 

,  El 

(152), 

EK 

( 152) s  A 

(152), 

2 

V 

(152), 

G 

(152) 

t  D 

(152), 

C 

(152) ,  X2 

(152), 

3 

X3 

(152), 

X4 

(152) 

,  X5 

(152), 

X6 

(152)  ,  X7 

(152), 

4 

SMI.A 

(152), 

SML8 

U52) 

,  SMLC 

(152), 

SMLD 

i 152) ,  SMLE 

(152), 

5 

EC 

( 152 ) » 

ER 

(152) 

•  SMLQ 

(152), 

SMLH 

(152 ; -  RIGA 

(152), 

6 

BIGB 

(152), 

CV 

(152) 

,  BC 

(152), 

BR 

(152) ♦  CHIC 

(152), 

7 

CHIR 

(152), 

CAPAC  ( 152 ) 

,  CAPAR 

(152), 

CRTC 

(152),  CRTR 

(152), 

8 

CRTPC 

(152), 

GOFR 

(152) 

,  FEW 

(152)  , 

CAR 

(152),  OKLM 

(  37) 

COMMON 

TELM 

(  37) 

,  EKLM 

(  37), 

ELM 

(  371,  FCLM 

(  37), 

1 

FRLM 

<  37), 

WLM 

(  37) 

,  QLM 

t  37), 

AMASNOi  371,  CHRNO  (  37), 

2 

ZP1 

(  37), 

ZP2 

(  37) 

,  SOLID 

(  37), 

ECHCK 

(  37),  RK 

(104), 

3 

RL 

(  37), 

RHOK 

(104) 

,  RDX 

(104), 

THETAKU04),  TEMP 

(  16), 

4 

HEAD 

(  12), 

MAXL 

,  MAXLM 

c* 

** 

C**************************4**************************Jjt***************** 
DIMENSION  CSQDU),  FM(l) ,  H4(l> ,  Y2U>»  XU)  t  YU),  HU),  SUMX2U) 
DIMENSION  SUMX3U),  SUMX4U),  SUHRHOU),  Q37<1),  Q38U),  H3U) 
DIMENSION  THE  TAT  (  23)  ,  KAMAX 123),  KCNTI23),  BUFS<750,1),  IDU2) 
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COMMON  /UNDLY/  HNU,SGNL, IHNU,NHNU,HNUP,NT, IM,IN,DHMU*THICK,NY 
C 

EQUIVALENCE  ISMLC,CSQQ1,  (SML8«EM)V  <SMLA,H4I,  <SMLO,Y2),  <GOFR,Xl 
EQU1  VALENCE  1 81 GA,  Y )  ,  C8 IG3 ,H1  ,  l  CRTR,  SUMX21  •  f  CHIC,  SUMX37  ,  (  8C  ,SUHX4) 
EQUIVALENCE  ISMLH,SUHRHOJ, (CAR,Q37J ,ICHIR,Q381,IBR,H3J 
EQUIVALENCE  ( ECHCK1 18 1 ,I0< 1)1 
C 

CAVEAT.  KAPPA  WILL  FAIL  ON  ANY  TWO  LAMBDA  REGIONS  OF  THE  SAME  MATERIAL 
C-  IF  THERE  IS  BETWEEN  THEM  AT  LEAST  ONE  REGION  OF  A  DIFFERENT  MATERI 

C  FOR  WHICH  A  DIANE  TAPE  IS  TO  BE  READ . 

C  ‘  AMASNOC 17-291  RESERVED  FOR  DESIGNATION  OF  DIANE  TAPE  UNIT. 

-100  FORMAT  1112,  12A6) 

IN=N 
IM*M  - 
NT  ■=  0 

,  NMAT  *  0 
11000  600  J*1,MAXLM 
IK*LMDAIJi 
IL*LM0A<J+1J-1 
IF  ( IL— INI  600,130,120 
220  IF  <IK-IN)  130,140,140 
1S(MK=IN 
-  GO  TO  170 

140  IF  ( IK— TM)  150,160,600 
150  IF  IIL-IMI  170,170,160 
140  IL«IM 
170  CONTINUE 

SOLI 01 201= IK 
SOLID! 21 J=IL 
SOLID! 22)* J 
L=GKLM( J1 

180  IF  {AMASNG<L*l?i I  190,190,320 
C  '  DONT  READ  A  DIANE  TAPE.  CALL  A  KAP  INSTEAD. 

190  GO  TO  ( 200 ,210,220, 230 ,240,250,260,270,280, 290 , 300 , 3 1 0 1  ,  L 
200  CALL  KAP1 
GO  TO  590 
220  CALL  KAP2 
GO  TO  590 
220  CALL  KAP3 
GO  TO  590 
230  CALL  KAP4 
,  GO  TO  590 
240  CALL  KAP5 
GO  TO  590 
250  CALL  KAP6 
GO  TO  590 
2 60  CALL  KAP7 
GO  TC  590 
270  CALL  KAP8 
GO  TO  590 


\ 
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290  CALL  KAP9 
60  TO  590 
290  CALL  XAPiO 
SO  TO  590 
300  CALL  KAPil 
SO  TO  590 
310  CALL  KAP12 
60  TO  590 

320  I  *  AKASN0(L+1?> 

IF  {NT— I I  330*460*330 
C  HEAD  DIANE  TAPE 
330  NT  «  ANASN0(L+17) 

NMAT  *  NNAT  ♦  1 

C  LIMIT  ON  NMAT  SET  BY  DIMENSION  OF  BUFS 
IF  INMAT-U  345.345.340 
340  SI  «  15.0931 
CALL  UNCLE 

345  IF  (KM AX I  370*350*370 
350  IF  (S15.EQ.C.1  GO  TO  420 

READ  (NT. 1001  IDA, ( ID( I ) .1*1*121 

REA0INT1  THNU.FNWO 

NMOxFNMD 

NHNU*THNU 

00  360  I  !*i*NHNU 

REAO(NY) 

360  REAO(NT) 

READ (NT)  HNU.FJAMAX 
4 AHA X*F 4AM AX 
IF( HNU.GT.O. )  GO  TO  365 
REAOINTI  (BUFSI I*NMAT).I«1*NH0) 

NHNU*1. 

REMIND  NT 
GO  TO  420 
365  Sl»15.0365 
CALL  UNCLE 

370  IF  (IHNU-ll  360.380.390 

380  READ  (NT.  ICO  I  IDA,(  10(11.1*1.  12) 

READ  (NT)  THNU.  FNWO 
NHNU  *  THNU 
NHO  *  FNMD 

390  READ  (NT)  HNU.  F4ANAX 
4  AM  AX  «  F4AHAX 

READ  f  NT)  (BUFSU.NMATI*  I»i.NWOI 
IF  ( IHNU-NHNU-ll  420,410,400 
400  SI  *  15.0957 
CALL  UNCLE 
410  REMIND  NT 
420  90  430  4A*i,23 
430  KCNTIJAi  *  S 
NCNT  *  1 
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00  450  JA*1,JAMAX 
THETATC JAJ  *  BUFS!  NCNT »NHAT  1 
KAMAXIJA!  *  BUFSINCNT+i.NMATl 
IF  ( JA— JAMAX)  440*450*450 
440  KCNT( JA+1)  *  KCNTIJAT  ♦  3  *  KAMAXIJA!  ♦  2 
450  NCNT  *  NCNT  ♦  3  *  KAMAXIJA!  ♦  2 
460  IL  *  SOLID! 21) 

IK  =  SOLID! 20) 

00  580  l=IK, IL 
TTST  *  Q374I! 

RTST  *  Q38! 1 2 
DO  470  JA*2* JAMAX 
IF  ITTST-THETATIJA!!  480*480*470 
470  CONTINUE 

JAT  *  JAMAX 
GO  TO  490 
480  JAT  *  JA 
490  KAM  *  KAMAXIJAT! 

KAMI  -  KAMAXI JAT-11 
DO  500  KA*2*KAM 
NCNT  =  KCNTIJAT!  ♦  3  *  KA  —  1 
IF  IRTST-BUFS1NCNT*3,NMAT3)  510*510,500 
500  CONTINUE 

KAT2  *  KA  -  1 
GO  TO  520 
510  KAT2  =  KA 
520  DO  530  KA=2*KAMi 

NCNT2  =  KCNTIJAT— 11  ♦  3  *  KA  -  1 
IF  (RTST-8UFS!NCNT2*3,NMAT ) }  540*540*530 
530  CONTINUE 

KAT1  *  KA  -  1 
GO  TO  550 
540  KAT1  *  KA 

550  NCNT1  *  KCNTIJAT-il  ♦  3  *  KAT1  -  1 
NCNT4  »  KCNTIJAT!  ♦  3  *  KAT2  -  l 
AA1  =  BUFS ! NCNT 1-2, NMAT) 

AA2  *  BUFSINCNTl+l.NMAT) 

AA3  *  BUFS INCNT4-2, NMAT,* 

AA4  *  BUFSINCNT4*1*NHATI 
AA5  =  BUFS ( NCNT 1* NMAT) 

AA6  =  BUFS! NCNT l+3*NMAT I 
AA7  =  8UFSINCNT4,NMAT! 

AA8  =  BUFS! NCNT4+3*NMAT! 

AA9  =  THETATIJAT) 

AA10  =  THETATIJAT-i) 

AA11  =  BUFS! NCNT 1-1* NMAT ) 

AA12  *  BUFS! NCNT1+2»NMAT I 
AA13  =  BUFSINCNT4-i*NMAT ) 

AAJ.4  =  BUFS!  NCNT4+2*NMAT) 

IF  ! AMIN1 ( A8S! AA6-AA5 ! »ABS( AA8-AA7 ) * A3S! AA10-AA9) ) l  570*560*570 
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560  Si  =  15.0915 
CALL  UNCLE 

570  AAC  *  AAI  ♦  (K7ST  -  AA5J  ♦  IAA2  —  AA1)  /  fAA6  -  AA5) 

AAE  =  AAI 1  ♦  (RTST  -  AA5)  *  4AA12  -  AA11I  /  ( AA6  -  AA5J 
AAO  =  AA3  RTST  -  AA7I  *  4AA4  -  AA3I  /  ( AAS  -  AA7J 
AAF  =  AA13  *i  RTST  -  AA7)  *  <AA14  -  AA131  /  i AAS  -  AA7I 

CAPARt  J)  =  EXP<  AAC  ♦  £ AAC  -  AACI  *  4TTST  -  AA101  /  4AA9  -  AA101I 

CAPACUI  =  EXPIAAE  *  £  AAF  -  AAE)  *  ITTST  -  AA10J  /  f  AA9  -  AA10H 

580  CONTINUE 
590  CONTINUE 
600  CONTINUE 
S15~6.0 

C  1NCUCE0  EMISSION  CORRECTION  FOR  KAPPAS  OBTAINEO  FROM  KA P  ROUTINES* 
C  NOT  DIANNE  TAPES,  IN  MULTI FREQUENCY  CASE 
IF  iKMAX)  610,650*610 
610  IFIHNUJ  650,650,615 
615  00  640  J=1,MAXLM 
IK  =  LHOA(J) 

IL  *  LHDAI J+l)  -  1 
L  =  OKLMUI 

IF  IAMASNOCL+17))  620,620,640 
620  DC  630  I~IK, IL 

630  CAPARl  I)=CAPARm*(  i.-EXPI -HNU/THETAII  i  M 

640  CONTINUE 

650  CALL  OVCHKiKOMYi 

GO  TG  1660,670),  KDMY 
660  SI  =  15.0152 
CALL  UNCLE 
670  RETURN 
END 
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Appendix  E 
PLNKUT 
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305 


VtUr*«t 


P 


i 

1  c 

!  SIBFTC  PLNKUT  FULIST,D£CK,RcF 

'  FUNCTION  PLKKUTI 8ETA1  *  SETA 2i 

s  C  THIS  ROUTINE  WAS  RECOMPILED  NOVEMBER  IS*  1964 

1  -  C 

DIMENSION  CSCD(l),  FMC II »  H4C1I,  Y2(ll,  Xlli,  YUJ,  HID,  5UMX2I II 
DIMENSION  SUMX3U),  SUMX4(1I,  SUMRHGIII,  Q37U),  Q38(l),  H3C 11 


c 

♦ 

c* 

S  P 

U  T  T  E 

R  CO 

MHO 

H 

** 

1  c* 

** 

i  c* 

** 

i 

* 

COMMON 

LMDAI 371# 

NR 

KSMLR  , 

IA 

9 

IB 

9 

ICA 

ICB  , 

1 

1 

KMAX 

9 

BLANXI, 

BLANK 2, 

BLANK.3, 

IAPI 

9 

IBP1 

9 

ICAPI 

ICBP1  , 

| 

2 

:  I 

9 

IG  • 

N RAD  • 

BLANK4, 

I  AMI 

9 

IBM1 

9 

2CAMI 

IC8MI  * 

\ 

3 

IIPl 

9 

IGM1  , 

I ALPHA, 

BLANK5, 

TH 

9 

TMAX 

9 

BLANK6 

OELPRT, 

\ 

4 

FREQ 

f 

CNTHAX, 

AR  , 

AS MLR  , 

PUSHA 

» 

PUSHB 

9 

60ILA 

BOILB  , 

i 

5 

CVA 

9 

CVB  , 

SLUG  , 

ALPHA  , 

HVA 

r 

HV8 

9 

HCA 

hce  , 

\ 

6 

EMINA 

f 

EMINB  . 

CA  , 

CB  , 

GA 

f 

GB 

9 

GL 

GR  , 

r 

t 

7 

RhOL 

9 

RHOR  , 

EPIO  , 

EPS*  , 

RIA 

f 

RZ8 

9 

RDIA 

RDIB  , 

8 

RPIA 

9 

RPIB  , 

RPDIA  , 

RPOIB  , 

TPRINT, 

TA 

9 

TB 

TC 

{  COMMON  TD  *  TE  ,  DTH2  ,  OTH2P  *  DTH1  ,  DTRMIN  *  DTK AX  * 

!  1  OTMAXl,  OTMAX2 ,  DTMAX3,  DTR  ,  SWITCH,  CO  , .CHIN  ,  DELTA  , 

1  2  GAMA  ,  WCRIT  ,  SIGMAQ,  AC  ,  AC03T4,  CNVRT  ,  SUMRA  ,  SUMRB  , 

\  n  n  t  »  nrt  t  i  m  nr\f  n  r.  n  r  n  n  1  rilP  r  t  r  ^  r* 


3 

ROIA 

,  RCIAM1, 

R0I8  ,  R0I8P1, 

GMS 

,  SI 

,  S2 

* 

S3 

4 

$4 

,  S5 

» 

S6  ,  S7  , 

S8 

,  S9 

,  SIC 

9 

Sll 

5 

S12 

,  S13 

» 

S14  ,  S15  , 

S16 

,  S17 

,  SIB  , 

S19 

6 

S20 

»  EO 

, 

FO  ,  TAU  , 

2ERO 

T  R 

(152), 

DEL  TAR ( 152) 

7 

ASQ 

1152), 

RO 

(152),  VD 

(152)  , 

RDO 

(152), 

SMLR 

(152) 

8 

DELR 

(  37), 

P 

(152),  PI 

(152), 

P8 

(152) , 

PB1 

(  152) 

COMMON 

P2 

(152),  SV 

(152)  , 

RHO 

(152) , 

THETA 

(152) 

I 

W 

(152), 

E 

(152),  El 

(152)  , 

EK 

(152) , 

A 

(152) 

2 

V 

(152), 

G 

( 152)  ,  D 

(152)  , 

C 

(152) , 

X2 

(  1523 

3 

X3 

1152), 

X4 

(152),  X5 

(152)  , 

X6 

(152) , 

X7 

(152) 

4 

SMLA 

(152), 

SML8 

(1523,  SMLC 

(152)  , 

SMLD 

(152) , 

SMLE 

(  152) 

5 

EC 

1152), 

ER 

(1523,  SMLQ 

(152)  , 

SMLH 

(152) , 

3IG  A 

(152) 

6 

8IGB 

(152), 

CV 

(152),  BC 

(152)  » 

8R 

(152) , 

CHIC 

( 152) 

7 

CHIR 

1152), 

CAPAC  ( 152 ) ,  CAPAR 

(152)  , 

CRTC 

(152), 

CRTR 

( 152) 

8 

CRTPC 

1152), 

GOFR 

( 152) »  FEW 

(152)  , 

CAR 

41521 , 

OKLM 

(  37) 

COMMON 

TELM 

(  37),  EKLM 

i  37), 

ELM 

(  37)  , 

FCLM 

(  37) 

I 

FRLM 

(  37), 

WLM 

(  37),  Qt M 

1  37  3  , 

AMASNO(  37), 

CHRNG 

l  37) 

2 

IP  1 

i  37), 

IPZ 

(  371,  SOLID 

(  37), 

ECHCK 

(  37), 

RX 

<  104) 

3 

RL 

{  37), 

RHQK 

( 104) ,  ROK 

( 1041 9 

THE! AK ( 104) * 

TEMP 

(  16) 

4 

HEAD 

<  12), 

MAXL 

,  MAXIM 

C 

c*  ** 

#******£*$***#$$#$***$$*«*#«***$***&*  *«£#*&^*»* ¥****♦*♦* 
EQUIVALENCE  (SMLCfCSGQ),  {  $MLB*FM)  ,  I  $MLA»H4 1  ,  (SMLD,Y2),  (GOFR.X) 
EQUIVALENCE  !6 IGA, V I » 1 8 1 GB ,H) , I CRTR , SUMX2I , i CHS C, SUMX3 ) ,  (6C,SUMX4) 


ISSSTItp'S  HAS  MAKC  1MftSF®g  “&AS  Fill®) 


EQUIVALENCE  t SHLH, SUWRHO ) , iCAR,Q37> , (CHiR,Q38J , { BR.H3 i 


C 

COMMON  /T8LC/  TBL(I2G1 


C 


C  WALKING/ //DO  NOT  COMPILE  AND  GO/// TABLES  MUST  8E  * 

C  SUPPLIED  * 


C»**»*»**»**** «***♦***♦♦♦****» *♦**♦»♦************»♦♦»**♦♦*♦»♦♦♦  »»♦»♦»»♦♦ 

C  ASSUME  BETA2  MORE  THAN  8ETAI 
IF  (BETA2  -  TSLIimO,  15,  15 

10  PLHKUT  -  .05132911  *  (BETA2**3  *'il.O-  8ETA2  *  1.375  -  SETA2  * 

1  .05)1  -  8£TA1**3  *  I 1.0-  3ETA1  *  1.375  -  SETA1  *  .051 1J 
REFURN 


15  J  =  1 

8ETA  =  BETA1 

17  IF  (BETA  -  TBLIllJ 18*  21.  21 

18  FCN  =  1.0-  .051329911  *  8ETA«3  *  (1.0-  BETA  *  1.375  -  BETA  * 

1  .0511 

GO  TO  35 

21  IF  (BETA  -  10.)  19,  30,  30 

19  I  =  1 

16  IF  ( 8ETA-TBLC 1+11 J  20,20,25 

20  FCN  =T8L( I+6C J+i SETft-TSL (I 1 1 *CT8U I+61J— TBL( 1+601 J 

1  /( T8L ( 1+1 1-TBL ( I 1 1 

GO  TG  35 
25  I  -  I  +  1 

IF  £  I  -  601  16,  30,  30 

30  FCN  =  .153989733  *  EXP(-BETA1  *  (((BETA  +  3.01  *  BETA  ♦  6.01  * 
1  BETA  +  6.0 
35  GO  TC  (37,  401,  J 
37  J  "  2 

OFCN  =  FCN 
BETA  =  BETA2 
GO  TC  17 

40  PLNKUT  =  OFCN  -  FCN 
RETURN 
END 
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SIBMAP 

♦ 

TBLC 

TBL 


PLNKTT 


PLNKTT  LIST  »OECK,REF 

THIS  ROUTINE  NAS  ASSEMBLED  NOVEMBER  9,  1964 
CONTRL  T8L,TBL+112 

DEC  1.0,1.04,1.08,1.12,1.16,1.20,1.26,1.32,1.38,1.44,1.50 

DEC  1.6, 1.7, 1.8, 1.9, 2.0, 2. 1,2. 2, 2.4, 2. 6, 2.9,3. 2,3. 5, 3.8.4.0 

OEC  4. 2 ,4. 4,4.6, 4. &, 5.0, 5. 2, 5.4 ,5. 6, 5. 8, 6.0, 6. 16, 6.32, 6. 48 

DEC  6. 64,6.80, 6. 96, 7. 12, 7. 28, 7. 44,7.60 ;7. 76,7. 92, 8.08, 8. 24 

DEC  8.40, 8. 56, 8. 72, 8. 88, 9.04, 9. 20, 9.36,9. 52, 9. 68, 9. 84, 10.0 

DEC  .965382,. 961696, .957807, .953716,. 949426, .944936, .937836 

DEC  . 930303, . 922348 , . 91 3982 • . 9C522C , . 88978 , . 87334, . 85598 

DEC  .83779, .81886,. 79926,. 77909, .73736,. 69435*. 62882,. 56372 

OEC  . 50062 , .44071 , . 40298 ,. 36721 ,. 33353 ,. 30201 , . 27266 , . 24547 

DEC  . 22041 , . 19740 ,. 17637, . 15722 , . 13984, . 12714, . 11544, . 10469 

OEC  .05482,. 08579,. 07753,. 06999,. 06311, .05686,. 05117,. 04601 

OEC  .04134, .03710, .03327, .02981 , .02669,. 02388, .02135 , .01907 

OEC  .01702, .01518,. 01354,. 01206,. 01074,. 00955 

ENTRY  PLNKTT 

ess  o 

END 
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